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1 ABSTRACT 

We suggest using two head-on magnetohydrodynamic (MHD) shocks to achieve 
supercritical nuclear fission in an axially elongated cylinder filled with UF 4 gas as an 
energy source for deep space missions. The motivation for each aspect of the design is 
explained and supported by theory and numerical simulations in this report. A subsequent 
report will provide detail on relevant experimental work to validate the concept. 

In this report we focus on the theory of and simulations for the proposed gas core reactor 
conceptual design from the onset of shock generations to the supercritical state achieved 
when the shocks collide. To model the shock generation and propagation, an MHD code 
developed by one of the authors at ISR is primarily used. The MHD model is coupled to a 
standard nuclear code (MCNP), primarily by the authors at the University of Florida, 
Institute for Nuclear Space Power and Propulsion Institute, to observe the neutron flux and 
fission power attributed to the supercritical state brought about by the shock collisions. 
Throughout the modeling, realistic parameters are used for the initial ambient gaseous state 
and currents to ensure a resulting super-critical state upon shock collisions. 

In some trials, power densities of 10 10 MeV.cm are seen at the time of maximum density, 
and this can grow to 1 0 MeV.cm or even greater if the gas is left unperturbed in the 
center of the shock tube. This is a novel effect that was not predicted and has not been fully 
evaluated as a possible means of fission power conversion; the cycle time should be on the 
order of a few milliseconds. 



NUMERICAL SIMULATION COUPLED MHD-MCNP CODES 

The broad aim of this research is to present methodologies and simulation results from 
research conducted into the modeling of shock wave driven gas core reactor systems. The 
initial concept under investigation could be described as a pulsed magnetic induction gas 
core reactor (PMI-GCR). The idea is that a strong ionizing set of twin shock waves would 
be generated in a fissile fuel filled shock tube, and that these shocks would then collide in 
the vicinity of a neutron moderator material providing a brief explosive burst of fission 
energy further ionizing the gas fuel and creating outgoing shock waves strong enough to 
allow an effect magneto-cumulative flux-compression generator (MCFG) system in pick- 
up solenoid coils, thus extracting pulsed electrical power. This energy extraction 
mechanism is the origin of the “PMI” in PMI-GCR. The reactor side of this system 
however is the primary component of interest in this report. Thus although the entire 
system involves three coupled stages, (1) shock generation, (2) shock interaction and 
fission energy release, (3) MCFG power extraction, the scope of this report extends only to 
the fission energy release stage with minor comments on the shock generation stage. The 
specific aim of this report is to document various simulation approaches to modeling the 
generation of the colliding shock waves and to present preliminary results showing the 
fission power energy generated as the two incoming shock waves collide. Power 
conversion (item (3)) is not included in this report, although previous reports have 
discussed the principle including Phase-I reports on this project. 

Unlike most pulsed high magnetic field research the present concept under investigation 
does not use destructive explosives and is intended to by a fully functioning power cycle 
design. Although the release of fission energy in a short burst, sufficient for net power 
production, seems likely to require destructive energies this need not be the case as the 
results of this report show. Indeed, with the disk-like geometry of two colliding gas shocks 
it is rather astounding that fission criticality can be achieved at all. The system is therefore 
quite safe, at least in that it requires a great deal of effort to produce supercritical 
conditions. The goal however is to demonstrate extreme shock-wave pile-up in order to 
demonstrate prompt criticality and useful energy release. To perfonn this feat in a power 
cycle is a tremendous challenge and may place demands upon materials that go beyond the 
current state of the art in shock tube technology and vessel integrity. The investigation of 
such materials issues and feasibility of cyclical operation of a putative PMI-GCR system is 
however beyond the scope of this report. This report investigates only the shock generation 
and fission power release on the assumption that materials will exist that can withstand 
many repeated firings of such a system. 
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This report summarizes recent work on developing modeling tools for simulating the 
pulsed magnetic induction (shock wave driven) gas core reactor concept (PMI-GCR). The 
goal of modeling the shock generation and fission power production for this type of shock- 
wave driven gas core reactor is a complicated task for many reasons. For one, there is no 
general-purpose code for coupled solutions of gas dynamic equations and neutron transport 
equations. Although such a code could be developed the software development time would 
be excessive. Earlier proposals sought ambitiously to model the entire system from shock 
generation to magnetic flux compression power extraction. Owing to difficulties in 
obtaining adequate general purpose nuclear and magneto-gas dynamic codes the current 
proposal is scaled back to just simulating shock generation and fission power production in 
the shock-tube reactor. Another difficulty concerns the highly nonlinear MHD shock 
evolution and interaction. For this reason a 1 -dimensional MHD code has been used and 
coupled to the Monte Carlo N-Particle (MCNP) code from LANL. In future work this will 
have to be replaced by 2D and eventually 3D MHD codes to estimate accurately the losses 
in the system, but for now the requirement of a working model to investigate various 
operational regimes and scaling principles for the PMI-GCR concept has required that 
these real effects be suppressed. 

This report then describes how the ID MHD code has been coupled with MCNP to yield a 
simulation tool suitable for exploring the ideal characteristics of a PMI-GCR system from 
shock generation initiation to fission energy release. The power extraction part of the cycle 
has not yet been implemented and requires solving difficult plasma kinetic problems as 
well as collection of preliminary experimental data on UF4 plasma characteristics before 
any realistic estimates of the system efficiency and specific power can be obtained. 

A Linux platfonn was chosen for compiling the codes, as it was felt that parallelization of 
the codes may be desired later. The available MHD code was written in Fortran, but 
input/output handling is easier to coordinate with C/C++ routines and Linux scripts. So a 
multiple code mentality was used, where separate executables would be run and 
coordinated by a script program. Although it is possible for a Fortran program to call a 
C/C++ program, and vise versa, it is easier to either write a single code so that data is 
accessible globally, or to write a script that runs separate executables. We choose the 
scripting approach. In this way the compilation of the executables is easier to manage 
because they are independent of each other. 
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This report also includes partial documentation for the code ‘mhdldshock’ used in ISR 
(Institute for Scientific Research (Fainnont, WV) and INSPI (Innovative Nuclear Space 
Power and Propulsion Institute, University of Florida) research into UF4-fuelled shock 
wave driven gas core nuclear reactors. The complete code uses a Linux script to batch 
process two main programs, ‘mhdldshock’ and ‘MCNP4C’. The first, ‘mhdldshock’ is a 
magneto-gas dynamic code that uses flux corrected transport (FCT) to solve the MHD 
governing equations, while MCNP is a Monte Carlo transport code used in the combined 
ISR-INSPI code to obtain nuclear data for input into the MHD code. The MHD code used 
is a “^’’-dimensional code which solves the initial-BVP equations for the variables 
(py x ,v v ,p,By). This is suitable for shock tube geometries where azimuthal symmetry holds; 
the axial magnetic field B x is constant. This enables qualitative scoping calculations for a 
shock wave driven gas core nuclear reactor design to be rapidly developed with a minimum 
of complication so that crucial experiments can be suggested. Once the physics of the 1- 
dimensional case are known and rough regimes determined for viable reactor operation 
(length scales, pressures, scalar gas conductivity, boundary conditions, and other dominant 
parameters) then full-blown 2 or 3D code can be deployed. This may be documented at a 
later stage of R&D, probably most feasibly only after some preliminary experimental data 
is gathered to support the preliminary modeling in this report. 

Sections 4, 5 and 6 are basically devoted to the theory of current driven shock waves and 
how they impact the conceptual design. For Section 6 of the report we review the 
conceptual design of a few possible PMI-GCR variations, where the PMI-GCR concept is 
regarded as a class in the wider category of ambient-sub critical pulsed reactor types. 
Section 7 of the report considers the task of coupling the MHD code to the general purpose 
Monte Carlo neutron transport code MCNP4C that was used to calculate fission power 
source terms and neutron flux in the shock tube and surrounding BeO moderator material. 
Subsections 7.1 to 7.4 describe the MHD model and deal with modifications to the MHD 
code to allow communication with MCNP. Subsection 7.7 describes the nuclear-MCNP 
modeling aspects of this project. 

After describing the coupling of the two codes simulation results are discussed in Section 8. 
Section 9 then discusses prospects for future research on the concept design for the PMI- 
GCR along with a final report summary for this phase of research on the power production 
leg of the system. 
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3 CURRENT DRIVEN SHOCK WAVES 

Several mechanisms were considered suitable for achieving localized ionization, 
acceleration as well as heating of the ambient gas in ways which minimize the damage to 
the boundaries, current driven shocks appeared most promising, because they fulfill the 
following conditions: 

3.1 Boundary Generated 

Shocks get triggered by the application of current at the boundary of the proposed reactor 
core. To be more precise, the so-called step function currents are used; these are localized 
currents with infinitely short rise times (or short compared with the characteristic time 
scales with a corresponding short spatial scale). Their highly localized nature will trigger 
localized gradients in hydrodynamic as well as electromagnetic parameters suitable in 
triggering shocks as well as exciting MHD fast wave modes, which can steepen to fonn 
shocks. 


3.2 Ionize the Ambient Gas 

The boundary-generated currents advect into the device along with other accompanying 
hydrodynamic variables. The electric field associated with the current can ionize the 
ambient UF4 gas as the shock travels into the medium. This is in addition to the shock 
heating ionization that can take place if the incoming wave steepens to suitably high 
amplitude. However localized ionization by currents will be most effective against high 
recombination rate in a problem of localized nature. 

3.3 End Wall Damage Minimized 

The use of currents as opposed to moving mechanical parts in generating shock waves 
insures minimum damage to the mechanical boundary; e.g., erosion, fatigue and interaction 
of nuclear material with the walls are potentially very serious issues for the boundary 
material, which should be minimized. 

3.4 Wave Steepening Capability 

Most hydro or even MHD shocks arise from discontinuities. As such in the shock frame 
steady state conditions apply, or in the lab frame a shock leaves a uniform state behind it. 
Such a unifonn state possesses global characteristics and will not be suitable for 
applications in which localized states with given hydrodynamic conditions are desired; e.g. 
conditions required for super-criticality. Only shock waves that arise from wave steepening 
can yield such characteristics, are time dependent, and therefore are not steady state 
solutions of the MHD or fluid equations. Current driven shocks can be used to fonn wave 
steepening. 


3.5 Why Two Shocks 

Two shocks are proposed as a means to enhance the steepening process. Consequently 
one could obtain from the two colliding shocks, one state that has density, pressure not a 
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linear sum of the incoming shocks. The advantage of this process is that one can achieve 
much higher densities and pressures at considerably smaller shock tube lengths. This will 
be most useful when dealing with applications that have serious constraints for the length 
of the device one intends to use. 


4 SHOCK GENERATION MECHANISM 

The device we propose to use relies on the application of a step function radial current at 
the boundary of a cylindrical device. Fundamentally we propose a design similar to that 
suggested by Patrick ( 1 959) 1 , Cloupeau (1963) 2 , Kholev and Poltavchenko (I960) 3 for 
shock studies, and more recently used for another application in plasma focusing by 
Moreno et al (2003) 4 . The device relies on generating a radial step function current by 
capacitor discharge. Depending on whether one uses a pulsed versus a steady boundary 
current, one has the option of triggering either a discontinuity driven or wave steepened 
shock, or a combination of both. In the applications above, however, the authors primarily 
used pulsed currents. 

To generate two shocks, radial currents are applied at the two ends. Each radial current Jr is 
associated with an azimuthal magnetic field Baz. Thus they result in a force Jr 0 Baz on 
the fluid on each end. The currents can be applied so that the force can be directed into the 
box at each end. Fig. 1 represents a schematic representation of the current-driven double 
shocks. If there exists an axial component of the magnetic field, the current Jr and its 
accompanying Baz will get transported into the cylinder; this is the so-called switch on fast 
shock. 


1 Patrick, Physics of fluids, 2, 589, (1959). 

2 Cloupeau, Physics of Fluids, 6, 679, (1963). 

3 Kholev and Poltavchenko, DokladyAkad. Nauk SSSR, Vol. 131, No. 5, pp. 1060-1063, April, 
1960. 

4 Moreno, Casanova, Corea and Clausse, Plasma phys. Control Fusion 45 (2003). 
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Figure 1. Schematic representation of the double shock generation at the boundaries using 
radial currents with their accompanying azimuthal magnetic fields. 


5 ANALYTIC MODEL 

For an elongated device as shown above, one can safely assume major spatial variations to 
take place along the axial direction. However, given the vital role played by the radial 
current and its accompanying azimuthal magnetic field and their advection into the device 
with the shock, the model should account for the variations of Baz and its accompanying 
flow velocity. As such a 1 .5 dimensional MHD model is most appropriate. 

5.1 MHD Model Equations 

We designate the axial direction to be x and the azimuthal direction y. The relevant MHD 
equations of continuity, momentum, energy and Faraday’s law are, respectively: 

dp/ dt + d/ dx(pV x ) = 0 

8(pV x ) / dt + d / dx{pV x 2 +P + 1 / 2 {B 2 - B x 2 )] = 0 

d(pV y ) / dt + d / dx[pV y V x - B y B x ] = 0 

d / dt[ 1 / 2 p(V x 2 + V y 2 ) + P l( r - 1) + 1 / 2 (B 2 + B 2 )] + 

dldx[pV x {M2{V 2 + v;) + yP l(j - 1 )p) + V X B 2 - V y B x B y ] = 0 

dB v I dt + d I dx[V x B v - V y B x ] = 0 

dB x ldx = 0 


5.2 The Geometry of the Model 

Fig. 1 shows the general field and current configurations in the model; Fig. 2 depicts the 
geometry used in obtaining the analytic as well as numerical solutions to the 1.5 
dimensional MHD equations above. 
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Figure 2. The field and current geometry used in the model. These represent the applied 
currents and their accompanying fields at the two ends of the proposed device. 


Here the applied boundary currents ( Jz ) at the two ends of the box trigger shock waves 
into the conceptual device by their respective forces ( J, x B ) forces. 

5.3 Analytic Solutions 

A shock can get triggered for each of the applied currents above. The most important 
characteristic of such a shock is the presence of the transverse component of the magnetic 
field as it travels into the conceptual device. Such shocks are labeled as switch on shocks 
and have been observed in both laboratory and space physical plasmas. In this section we 
will derive the conditions under which such shocks can get triggered. Although ultimately 
the method used to generate the shocks in the simulations will be different, they are based 
on the switch on shock concept. The switch-on shock, though in its simplest form results 
from the application of a single step function current pulse, has exact analytic solution. 
That can provide useful insights into the simulation results, which will follow. 

The switch on shock solutions are obtained using discontinuity theory. Using the standard 
wave frame calculations, as in Kemp and Petschek 5 and Bazer 6 , one will obtain the 
following solutions to the above MHD equations after applying the current and field 
discontinuities above as follows: 


K 2 2 =B x 2 /p 2 =b Xi 2 

(PiK, /2)( r + i)/( r -i)F X 2 2 +[- r /(r-i)(p 1 K I 2 + B 1 -B y 2 !2)-B y 2 ]v Xi 
+ 1 / 2 Pl V Xi 3 +y/(y-l)P l V Xi +l/2B x 2 B yi 2 /( Pi V Xi ) = 0 
Vy 2 = B x /( Pl V Xi )B yi 


5 Kemp and Petschek, Physics of fluids, 2, 599, (1959). 

6 Bazer, Astrophys. J. 128, 686, (1958). 
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where the subscripts 1 and 2 refer to the un-shocked and shocked regions respectively. 
Clearly B = 0 , and B = B from the last equation above. Using V x and V y of the first 
and the third equations into the second equation gives the following relationship for the 
switch on field component: 

B yi 2 /B x 2 =2(p 2 /p l -l){l-.5(r-l)(p 2 /p l -i)-rP l /B x 2 } . 

This equation furnishes both the condition of the onset as well as some important 
implications of the switch on shock; i.e., for the right hand side to be positive we must 
have: 


2 2 

(a) p 2 t p x > 1 , which also implies that V / b > 1 ; 

this is consistent with steepening. Furthermore, since V > b (the Alfven speed in the 
upstream un-shocked region) then this shock is a fast shock. Since the flow velocity 
downstream V is simply Alfven speed, this shock must also be an intennediate shock. 

The second condition that results from the above is: 

(b) p 2 / A = V Xi 2 /b Xi 2 < (r + 1)(^1 2 / 2)[( 7 - \)M 2 / 2 + 1] ; 

where /V/, 2 = V x 2 KyP ] I p x ) is the Mach number in the un-shocked region. This 
relationship’s right hand side is simply p 2 / p x of a pure hydro shock. Therefore, as 
expected for the MHD shock, the density steepening is lower than its corresponding hydro 
shock because some of the shock energy goes into the magnetic field steepening. However, 
as will be seen in the simulation section for the double shock simulations, heating can be 
substantially higher in this case, compared to the hydro case, if the magnetic energy that 
goes into the MHD shock can be turned into thennal energy by shock collisions. 

The equation above for B v 2 / B x suggests shock generation by perturbation of only one 
parameter, B y . Most discontinuity generated shocks demand perturbation in all the 
independent parameters. This property has far reaching implications: it simplifies shock 
generation, it enables calculation of all the downstream parameters in terms of one, and can 
be used to trigger non steady state and time dependent downstream states. 

The B / B x equation suggests at least two means to manipulate downstream parameters: 
(a) by application of very large currents with the resulting downstream parameters of 
desired size, or (b) by holding a relatively small current steady for a finite length of time. 
Case (a) corresponds to a pulsed current case, while case (b), which will demand 
maintaining open boundaries while the current is on, resembles the mechanism of a thruster 
into the conceptual tube at each end. From this point on we will refer to the two cases as the 
pulsed versus magnetoplasma dynamics (MPD) thrusters respectively. 
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In the former case, the size of p 2 and consequently V and V are completely detennined 
by B v and since the boundary closes as the current is turned off, instead of a discontinuity 
with upstream and downstream states possessing global characteristics, a single pulse 
results (with front speed V ) followed by a rarefaction wave. In the latter case, the 
continuous current pulse can be thought to result from infinitely many short duration 
pulses. There each pulse enters a denser medium perturbed by its predecessor with larger 
speed. It can, therefore, take over and merge the predecessor and give rise to a larger pulse. 
In both cases the impact is time dependent, strongly governed by the steepening process. 

We shall next present the numerical simulation results of the MHD code as well as the 
MHD-MCNP coupled codes modeling a conceptual reactor from the shock generation to 
the fission triggering. 


6 PMI-GCR CONCEPTUAL DESIGN OVERVIEW 

Generally for the shock generation a coaxial arrangement of anode and cathode as depicted 
in Fig. 3 can be used. Such a design ensures very strong self field ( B e ) accompanying a 
strong radial current J and therefore strong pondromotive ( J x B e ) force. The strong radial 
current J serves to ionize the gas too. 

Where the radial current ends the field of the solenoid begins which can trigger the switch- 
on shock and allow advection of both J and B e into the device. Application of two such 
currents (with mirror symmetry) at the two ends of a coaxial cylindrical device will result 
in head-on shocks. 

Although many variations on the same theme could be generated, this report is concerned 
only with two specific classes of pulsed magnetic field shock-driven reactors. One is 
tenned an “MPD compressor mode” device; the other is a “pulsed high magnetic field” 
device. The distinction between the two types is the method of shock wave generation at 
each end of the shock tube reactor. In practice the distinction between these two shock 
generation methods (“compressor-type” and “pulsed-type”) becomes blurred and 
indistinguishable when the pulse time is set fairly large and the current density at the 
boundary is reduced: in that case, a pulse becomes almost indistinguishable from a short- 
duration MPD-compressor mode. The reason a sharp distinction cannot always be made is 
because, to operate in a pulsed high 5-field mode, the gas fuel has to be highly ionized, and 
this generally requires a large current discharge through the gas. Consequently 
ponderomotive force effects will play a dominant role in addition to the desired magnetic 
discontinuity shock inducing effect. Thus, for high 5-fields, shock waves will be created in 
the gas through two effects, one through the ponderomotive force and the other via 
magnetic and the resulting pressure discontinuity. These two effects may compete or 
cooperate depending upon the design of the shock generator, its geometry and electrode 
configuration. 
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The term “MPD compressor mode” is used because for smaller current densities and longer 
pulse duration times, the shock generation arises mainly from nonlinear wave steepening as 
described in the analytic section above and Wu ^ 7] . For short pulse times and large current 
densities the shock generation can be said to arise more from steep magnetic pressure 
gradient at the boundaries giving rise to shocks from the feedback of magnetic pressure 
discontinuity into the gas through classical magneto-gas dynamic coupling, as described in 
Knoepfel [S| . There is, however, a regime where a true pulsed mode could be identified as 
distinct from a compressor mode by examining the time scales for collisional, relaxation 
and ionization processes; such an analysis should reveal when the transition from 
discontinuity-induced to pondermotive force-induced shock occurs. This report did not 
look closely at these issues however; since we believe that both mechanisms are at work, 
but with pondromotive to dominate the fonner. In both cases the coaxial setup depicted can 
be used with simple change in the duration of the applied currents and UF4 gas injection. 



Figure 3. Representative generic coaxial or annular electromagnetic-driven shock tube. 

The inner electrode here can be extended the entire length of the tube if a steady state annular shock wave is 
desired. This is probably unsuitable for the PMI-GCR system, however if the electrode was extended only part- 
way into the shock tube an annular tube geometry could be maintained by making the central tube ut of some 
moderator material perhaps even mixed with a small fraction of solid fissile fuel . 


[ * * 7 ] C.C. Wu, “Formation, Structure, and Stability of MHD Intermediate Shocks”, J. Geophysical 

Research, 95, N0.A6, 8149-81475, 1990. 

[ 8 ] H.Knoepfel, “Pulsed High Magnetic Fields”, North-Holland, 1970, in particular chapters 6 to 

8 . 
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7 SHOCK GENERATION AND INTERACTION MODEL 

The overall viability of the PMI-GCR concept depends upon three crucial subcomponents: 
(1) the efficiency and efficacy of generating strong shock waves in the comparatively large 
diameter, high pressure shock tube, (2) design of the core interaction region for release of 
fission energy, and (3) magneto-cumulative flux-compression power extraction. This 
section discusses the first of these components. 

7.1 MHD Equations for Electromagnetic Shock Tube 

The governing equations of weak MHD can be derived by adding electromagnetic source 
terms to the Navier-Stokes equations. The continuity equation is unchanged because the 
mass density of the electromagnetic field is negligible. We can write mass conservation in 
the following fonn. 


^ + V-(pv) = 0 
dt 

Here p is the fluid density, and v the velocity, so pv is the momentum per unit volume of 
fluid. The momentum equation is a vector equation for pv. If we apply Newton's 2nd law 
to a control volume, with external forces per unit volume F ext , and pressure-stress tensor , 
and then apply Gauss’ law to reduce the integral equations to differential fonn, then 
momentum conservation takes the fonn, 

^(pv) + V-(pvv) = F 6x1 + V • T 
dt 

In this research the only body force is taken to be the ponderomotive force per unit mass, 
F e ' n , which is derived from the Lorenz force. The Lorenz force law on moving charges 
together with the charge conservation law imply that the external body force is therefore, 

F em =p {c) E + -JxB 
c 

here p Ui is the charge density, and E, B and J are the fields and total current density 
respectively in the stationary observer frame that sees the fluid moving at velocity v. For 
overall neutral gases one can safely assume p {c) = 0 over even small volume elements. The 
choice of units here is Gaussian (cgs), for which po=so~l. We also make the usual MHD 
approximations (scalar conductivity, negligible Hall term in Ohm’s law), for which 
Ohm’s law reduces to, 

J = oE + — vxB , 
c 

so that cr is the scalar electrical conductivity of the fluid. 
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An energy equation for MHD can be derived easily by assuming a linear isotropic fluid, 
and then balancing (right-hand side say=) the rate of increase of total energy with the sum 
of (l)work done by surface stresses, (2) work done by external body forces, (3) heat 
conducted into the system, (3) energy radiated into the system (= left-hand side). In this 
research heat conduction and viscous stresses are ignored. The energy equation can be 
formulated either by conceptualizing the electromagnetic field as part of the system in 
which case there is no external body force. However one must consider the electromagnetic 
field energy on the left of the energy balance and the flux of electromagnetic energy 
c !{An)V • (E x B) radiated away (therefore subtracted on the right hand side of the energy 
balance). Alternatively, one can conceptualize the MHD problem as two subsystems, a 
mechanical part and a field part, and derive an energy balance on the fluid only (however 
then one ignores field energy density and flux), but instead one must calculate the work 
done by the (now external!) ponderomotive force, c -1 v • J x B , as well as the Joule heating 
= j / <j . Strictly speaking one should use the conduction current density J' seen by the gas 
to calculate the Joule heating, but for our purposes J ~J. Both terms must be added to the 
right-hand side of the energy balance. 

Thus, there are two entirely distinct ways to conceive of conserving energy for the entire 
MHD system, but both turn out to give an equivalent energy equation. We prefer to write it 
using the unified system fonn, which looks at field energy density and electromagnetic 
flux. For reference, the full equation (though ignoring line radiation of the gas) would be, 

— (pe + U (em) ) + V • (pe\ - pi • v) = V • (t ( ™ c) • v) - — V • (E x B) + V • (xVr) . 

dt 4 n 

Here there are no additional forces doing work on the system: e is the specific total 
mechanical energy of the gas, e=u(T)+0.5v 2 , f/ fr '" ) =(E 2 +B 2 )/8rc~B 2 /87i is the internal 
energy density of the field, p is the gas pressure, x (mt) is the viscous stress tensor, and /vis 
the thermal conductivity. 

To form a closed system of equations for (p,p,\, B) an equation of state of some form is 
assumed to calculate the temperature form (p,p) and the internal energy of the gas is taken 
to be u""=U()+c v (T-To), where uo= u mt (To ) for some reference temperature. For UF 4 fuel 
M UF4 =0. 3 140225 kg.moF 3 , and y~1.07, so cv~/?ur 4 /(y- 1 )=378.247 J.kg '.K 1 = 3.78xlO () 
erg.g .KT 1 . With appropriate boundary conditions assumed, there are five equations above 
for eight unknowns. Three equations can be derived from Maxwell’s homogeneous 
equation for dB/dt. These are the B-field transport equations. Using the generalized 
Ohm’s law approximation, as well as using Maxwell’s equation for the current (and 
ignoring the negligibly small displacement current), so that 

J = — VxB, 

An 


one ends up with, 
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r)R 

— = Vx(vxB) 
dt 


Vx(VxB) 

4 n 

or because V x (V x a) = V(V • a) - V“a for any vector a, and noting V • B = 0 , 

— = Vx(vxB) + — V 2 B. 
dt 4 n 

This completes the MHD equations. For application to shock tube simulations these 
eight equations can be reduced to a set of five coupled PDEs for (p ,p,v x ,v y ,B y ). 

7.2 Electromagnetic Shock Tube System 

The 1D-MHD code that has been coupled to MCNP currently employs a current sheet 
(which could be either formed (1) radially in practice from coaxial electrodes, or (2) as a 
planar sheet between vertically opposite anode and cathode) and crossed transverse 
magnetic field at the boundaries of the neutron moderator region to launch MHD waves 
into the moderator region from both left and right boundaries of the shock tube. These 
waves build up in intensity until the current and transverse magnetic fields are turned off 
after some designated time. The details of this gasdynamic/MHD shock production were 
treated above in the analytic section. The shock production mechanism is however not 
completely irrelevant to the nuclear modeling. There may be reasons why a particular 
shock generator will not be suitable. For example, a diaphragm shock tube is too difficult 
to set up for continuous cycling, and a coaxial quasi- 1 -dimensional shock tube will not 
allow efficient fission power production, but a coaxial electrode generator may be used if it 
can generate a shock strong enough to propagate into a proper gas-filled shock tube region. 
Many of these issues give rise to questions that can only be answered by at least a 2D 
magnetogasdynamic simulation, and possibly only by supplementary experimental data. 
The final project report will eventually comment more upon these questions. 

To simulate the shock formation and interaction for a conceptual shock-driven gas core 
nuclear reactor, the governing single fluid MHD equations were reduced to cylindrically 
symmetric form with further simplifying assumptions amounting to setting the axial 
magnetic field f? r =constant. To reduce further to “1 ^’’-dimensional geometry Cartesian 
coordinates were used to set up the discharge current and associated driving magnetic 
intensity. Thus, to simulate either a quasi- 1 -dimensional coaxial shock generator or other 
Lorentz force type MPD compressor, one needs crossed magnetic field and current sheets, 
i.e., at the boundary one requires B perpendicular to J and that both are perpendicular to the 
desired axial flow direction of v r . Thus, choosing a convention the simulation had B z =v z = 0, 
and J y = 0. The semi-continuous discharge current J z is then set by using Maxwell’s 
equation and ignoring the negligibly small displacement current, so that 
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J = — VxB, 

An 

which is used to set whatever B=(0,Z? > ,0) at the tube boundary is required to obtain a 
given desired current J : (hence, ponderomotive force). 

For the shock tube symmetry then the derivatives dldy,dldz^> 0, f?_ — » 0 and 
/? A =con slant. Furthermore, in the simulations to be discussed in this report, negligible 
thermal conductivity and viscosity are assumed, and radiative losses are also not modeled. 
Under these assumptions the governing MHD equations become, 


dp d 


8 d 

t;O x ) + — 
ot dx 


— (pv v ) + — 
dt y dx 


pv 2 +P + P(B 2 y -B 2 ) 
8 n 


= 0 , 


P v y v x ~^( B y B x) 


= 0 , 


pe + ^-{B 2 +B;) 
an 


d_ 
dt 

8By C 

dt dx 


+ - 


dx 


1 


pv x (e + P/p) + —(v x B y - v y B x B y ) 


{An) 2 8x\ 


'By dB y 
a dx 


( v y B x -v x B y )+ 


- 2 d 2 B . 


Ana dx 2 


This form is fairly easily to put into discretized form. However, by converting to a 
rationalized Gaussian system of units the constants c and An can be eliminated from the 
equations thereby allowing an even more convenient discretization for computer coding. 
Specifically, looking at the energy equation it is clear that defining, 


b = 


B 

An ’ 


we can eliminate 4 ti. Then, by setting, 

An 

we can entirely eliminate the constant An from all the equations. Note that this amounts 
to also redefining the electric field as E, and current density as %, which can be computed 
directly from b and q using the Gaussian formulae with An — » 1 , and c — > 1 , and can be 
related to the conventional Gaussian system E and J by the following relations, 
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c 



E, 


X = 



J, 


For example, the ponderomotive force law, F em = c 1 J x B , (for a neutral plasma) then 
becomes just, F e ' n =yx b , in the new rationalized units. Note that in converting to 
these units the only drastic change is that the electrical conductivity s has different 
dimensions, in Gaussian units [cr]=s , whereas in the new system [q]=cm .s, or 
equivalently one can consider the resistivity which has dimensions of Time in the 

cgs system, now has dimensions of [Length - Time ] or units of cm s in the new system 
where c~ and 4n do not appear. As a quick example, consider the scalar electrical 
conductivity of a steady state UF 4 plasma at P= 1 atm, Z , =2000°K roughly, then cr~l .0 
[Siemens nT 1 ] from data in MKS units, which is cr «8.98xl0 9 [s _1 ] in Gaussian (cgs) 
units. But in the new rationalized system used for computer coding we would use q = 
8.98x10 ' 47 t/c = 1.69x10 " [cm - .s] in the new units. In other words the electrical 
resistivity that needs to be used for the computer code would need to be set about 
c 2 /4ti ~5.3x 1 0 14 times larger than the conventional electrical resistivity of the plasma in 
Gaussian units. With these modifications taken into account, the system of equations for 
computer discretization in conservative form are further simplified. 


The key thing to note is that the non-electrodynamic variables (x,t,p,\,e,P) are still in 
Gaussian units of g, cm, s. The new variables only affect the appearance and dimensions 
of the electrodynamic terms. As a final note before describing the computer code used to 
solve these equations, it should be remarked that strong shocks are capable of ionizing 
gases, and across electrical discharges in gases the conductivity approaches that of a fully 
ionized plasma: therefore, it is not too extreme, for the application that this report is 
concerned with, to make the approximation of infinite electrical conductivity, then 
Mg — » 0 : hence, the diffusive terms in the last two equations may be ignored as a rough 
approximation. Of course simulations hoping to capture energy efficiencies for the 
processes involved obviously cannot make such a rash approximation; indeed, such goals 
are not really compatible with neglect of thermal conductivity, viscosity and radiative 
losses. Such features of realistic shock tube processes are, however, deferred to possible 
later stages of research for the current project. 
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7.3 Solution of the MHD Shock Tube Equations 

Although the report does not cover the numerical method used to solve the above 
equations. However some brief comments will be given where appropriate to note the 
particularities of the numerical algorithm written in Fortran code that was used to obtain the 
results reported on in Section 8. One complication in particular frustrated efforts to use a 
full flux-corrected transport MHD solver algorithm. The MCNP code consumes 
considerable CPU time each time it is called, therefore to limit simulation run times to 
days rather than weeks the MHD code was allowed to run for a given number of time steps 
and/or simulation physics seconds before pausing to call MCNP for updating the fission 
power and neutron flux source terms. In addition, for adequate convergence in the MCNP 
tallies, the shock density profile had to be collapsed into just a few cells, where the shock 
waves gave steep density profiles the MCNP file input writer script made the collapsed 
cells finer, but where the density varied slowly many cells of the fine MHD grid mesh 
would be merged into one or two large cells. This made the fission power density q "' , and 
neutron flux ®„, calculations rather crude. It was decided therefore that the fission source 
tenn in the energy equation would not be added to the MHD solver algorithm until the 
coupled codes could be verified as accurately reproducing the MHD-shock simulation 
profiles for a non-fissioning gas. The intent was to proceed to include the fission power 
source tenn (as well as the fission-fragment and neutron flux ionization enhancements to 
the electrical conductivity) but all efforts to debug and validate the modified MHD code 
and script programs failed to culminate in time to allow a further round of debugging and 
testing of the fully coupled MHD-MCNP code with fission source term feedback. 


17 



A leapfrog trapezoidal algorithm (LFT) was used for transporting the fluid variables a 
given time step. The main difference between gasdynamic LFT and the MHD LFT used in 
this report is that the MHD system includes additional wave speeds, the fast and slow 
magnetoacoustic waves, as well as Alfven waves. Details of leapfrog trapezoidal 
algorithms for MHD and flux-corrected transport methods can be found in the references 
[ 9 , 10> n, 12^ 13 j -phere is also a time scale characteristic of magnetic field diffusion that 
affects the stability of the solution. The fast wave velocities along with the gas dynamic 
acoustic wave velocity were used to detennine the allowable maximum time step dt\ for a 
fixed spatial mesh size dx, consistent with CFL stability. The time scale for magnetic field 
diffusion was then set as dt 2 = 0 . 5 dx !q , where q is the electrical resistivity defined above. 
The time step for each iteration of the MHD equation solver was then taken as dt = min (dt 2 , 
dt\). 

The remainder of this section considers the task of coupling the MHD code to the general 
purpose Monte Carlo neutron transport code MCNP 4 C that was used to calculate fission 
power source terms and neutron flux in the shock tube and surrounding BeO moderator 
material. 


9[ 9 ] D.L. Book, J.P. Boris, and K. Hain,, “Flux-Corrected Transport. I. SHASTA, A Fluid 
transport Algorithm that Works”, J. Comp. Phys., 11, 38-69, 1973. 

[ 10 ] D.L. Book, J.P. Boris, and K. Hain,, “Flux-Corrected Transport. II. Generalizations of the 
Method”, J. Comp. Phys., 18, 248-283, 1975. 

f 11 ] J.P. Boris, D.L. Book,, “Flux-Corrected Transport. III. Minimal-Error FCT Algorithms”, J. 
Comp.Phys., 20,397-431,1976. 

[ 12 ] S.T. Zalesak, “Fully Multidimensional Flux-Corrected Transport Algorithms for Fluids”, J. 
Comp. Phys., 31 , 335-362, 1979. 

[ 13 ] C.R. DeVore, “Flux-Corrected Transport Techniques for Multidimensional Compressible 
Magnetohydrodynamics”, J. Comp. Phys., 92, 142-160, 1991. 
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7.4 Modifications to the MHD Code for Serialization with 
MCNP 

The general purpose Monte Carlo code MCNP4C is suited for high precision physical 
modeling of fundamental radiation interaction phenomena. It is less suited for dynamic 
system simulation, where a discrete ordinance code might be favored. Nevertheless, in the 
nuclear engineering community MCNP is a favored tool; hence, it was chosen for use in 
this project. The main difficulty in using MCNP for the shock-driven reactor system 
modeling is that it takes many minutes for MCNP tallies (Monte Carlo schemes all involve 
tracking individual particles and building up statistics from such primitive numerical trials) 
to converge to within reasonable tolerances. But the MHD simulation code updates the gas 
pressure, temperature and density roughly every microsecond. Yet the shock waves travel 
at finite velocity, hence for accurate fission power calculations the MCNP code needs to be 
called many times during a reactor simulation cycle. With some 2000 to 10,000 time steps 
typical for the Leapfrog trapezoidal MHD solver it immediately dawns upon the 
simulationist that it will be untenable to update the fission power source terms every step of 
the MHD shock evolution. Therefore, some criterion, or set of criteria, is needed to 
detennine judiciously when to pause the MHD code to call MCNP to update the fission 
power terms. For this report a combination of a minimum lower gas density (below which 
fission events could be regarded as negligible) and a minimum time in seconds over which 
the density profile would not change too rapidly and yet exceeding an average neutron 
lifetime in the gas were the chosen criteria. However, even these conditions often resulted 
in excessively slow simulations, so a third criterion was used, which was to explicitly limit 
the maximum number of times that MCNP is called to a fixed number over the entire 
simulation cycle. 
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Figure 4. Flow diagram of operating script to control simulation linking MHD1 DSHOCK and 
MCNP. 

This logical loop is maintained until the MHD1DSHOCK code completes the interaction of the two shock waves 
or reaches its maximum assigned integer time step (or computer memory is saturated). 


Serialization of the MHD and MCNP codes, therefore, involved writing a Linux script 
program to coordinate the following executable tasks: (1) starting the MHD code and 

(2) running it until the conditions required for calling MNCP are attained, then 

(3) preparing input cards for MCNP on-the-fly using intermediate program codes, then 

(4) running MCNP, (5) checking for convergence, if it hasn’t converged then (6.1) keep 
running MCNP, then (6.2) upon convergence preparing an input card for the MHD code to 
(7) resume from the previous time step, and repeating steps (2)-(7) until the MHD code 
finished a cycle (reached it’s maximum set time step.) A flow diagram of the script logic is 
shown in Figure 4. 
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7.5 Getting the Coupled Codes Started 

The most acute sensitivities of the simulation results to input parameters are the electrical 
conductivity and neutron source strength. Low conductivity, or high resistivity poses a 
problem because it sets the time scale for magnetic field diffusion. If the magnetic field is 
large at the boundary grid cells then sudden changes in intensity can lead to anomalous 
numerical results such as negative pressures being produced. The neutron source strength 
is adjustable and easy to deal with, so if either the fission power density or neutron flux 
exceeds the compiler single precision maximum for floating point numbers then (implying 
an unwanted highly supercritical core) one can lower the ambient source strength. These 
peculiar problems will be dealt with in Section 7.6. For now we will discuss the practical 
details of finding good starting conditions. 

As of the present date the codes function together well, but the user must exercise care in 
selecting the input parameters described above in Appendix Table 7. One must ensure that 
the MCNP pre-processor does not collapse the data from the MHD fine grid mesh into a 
single cell, because this can cause problems for the MCNP post-processor or MHD pre- 
processor. This can happen if MCNP is called early on in a simulation before the gas 
density starts to get perturbed by the shock generator. This problem gets amplified when 
the MHD simulation time step is very small and of the order of 10’s of nanoseconds. With 
such very small time steps the shock generator has not had time to produce significant 
shock wave build-up of gas density, so the pre-processor is in danger of collapsing the grid 
to a single cell. This will be amended in future versions of the scripts. However it is 
nevertheless interesting to examine how such small time steps can arise, because they are 
symptomatic of another problem: resistive shock instabilities. That is, the resistive shocks 
at high Mach numbers encounter physical instabilities. In the actual physical problem, as a 
wave steepens and its scale becomes smaller than the resistive scale length, other 
dissipations with shorter scale lengths will be needed to prevent the nonlinear growth of the 
wave. So either viscosity or other transport will be needed or in doing numerical 
simulations numerical diffusion can be increased by increasing the grid size or other means. 

Section 7.6 will discuss some of these peculiar problems. General comments on how to get 
the codes started and running together will now be summarized. 
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7.5.1 Finding Sensible Initial Conditions and Inputs 

Physical sensibility and research interests should determine ambient pressure and density 
“Pin”, and “rhoin” in Appendix Table 7. Magnetic fields cannot exceed megagauss 
intensity (and even then only for short pulsed fields) so there are physical limits that 
determine reasonable shock-drive boundary conditions “byi” in Appendix Table 7. The 
neutron source strength cannot reasonably be much in excess of 10 2<) neutrons per sec, 
because even accelerator driven systems would be pushing practical limits to provide such 
and enonnous quantities of neutrons. It takes anywhere between 5 to 20 minutes to run one 
call to MCNP, so is unreasonable to set “ncalhncnp” to much more than 100 to 200 calls 
for rapid benchmark research. Given (as discussed later) the number of time steps required 
to see the two shocks collide can be in the 10’s of thousands, one clearly should limit the 
total number of times MCNP will be called by basing calls to MCNP upon either real 
(simulation) time in seconds, or upon ncalhncnp. One can always run the coupled code 
script and use the system’s keyboard kill command to quit a simulation early if it looks like 
the shock collision has been captured. 

Once a desired et of input parameters has been selected the best way to run a simulation is 
to first turn MCNP off by setting “ismcnp=F” in the namelist file Tnhdini.dat’. One can set 
a very large maximum time step, say “nend= 100000” in this case, and wait a few minutes 
to see the results. Note that with MCNP thus turned off it is simpler to run the MHD code 
without going through the Linux shell script, just run the native MHD executable. 

If negative pressures which are a symptom of resistive shock instabilities, are encountered 
then, generally, if one wants to keep the same shock tube length, resistivity and grid 
resolution, then one needs to either, (i) artificially reduce the maximum allowable time step 
by reducing the fraction “diff_tscale”, or, (ii) decrease the applied magnetic field by 
decreasing “byi”, to reduce the magnetic Reynolds number. It’s a matter of trial and error. 
However, by choosing option (i) to decrease diff tscale, one may in effect reduce the time 
step dt set by the subroutine SETDT to such a small value (maybe even nanoseconds) that 
the simulation stops before the shocks collide. To avoid missing the shock collision one 
has to inspect the simulation (again with MCNP turned off) to see if with the smaller time 
steps the simulation does indeed advance far enough to produce a collision, if not then one 
must increase “nend” to get more total time steps. 

If no negative pressures are seen (and of course of sufficiently strong shocks were 
produced) then the inputs are probably fine and one can set “nend” to roughly twice the 
time step at which one saw the shocks collide (to get a decent history of the post collision 
evolution). Then turn MCNP on by setting “ismcnp=T” and run the shell script again and 
wait for the full results of the coupled MHD-MCNP codes. 

For smoother (less discrete) fission power and neutron flux data one can increase the 
number of times MCNP is called by either increasing “ncalhncnp” or decreasing “dtmcnp”. 
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7.6 Trouble shooting the Simulation Codes 

The peculiar problems encountered while testing the MHD1DSHOCK and MCNP4C 
coupling are now described. 

7.6.1 MCNP4C Convergence Criteria 

To force MCNP to converge rapidly one can do two things, 

(1.1) One could collapse a lot of cells from the MHD mesh to just a few, averaging out 
the gas density. This will mean MCNP will not have to wait too long for sufficient tallies to 
build up in each cell. To do this edit the file ’mcnpprob’ and change the ninth number in the 
line, 


mcnpspec 1 1 20 0 0.005 0.05 0.05 40 1.0 62.5 

the fifth number here is the maximum permissible standard deviation for k e f that MCNP 
uses to detennine convergence. Ideally for high statistical confidence, the maximum 
standard deviation should really be around 0.005 or less, however smaller values will 
drastically slow down the simulation (by requiring more cycles be simulated in MCNP 
for better statistics). The sixth and seventh numbers are the maximum permissible 
standard deviation in the flux and fission heating tallies respectively for individual cells. 

(1.2) One can also specify a statistically large maximum permissible standard deviation for 
k e ff. This is done by editing the file ’mcnpprob’. At the line, 

mcnpspec 1 1 20 0 0.005 0.05 0.05 40 1.0 62.5 

the fifth number here is the minimum standard deviation for fy// that MCNP uses to 
determine convergence. Ideally for high statistical confidence the min. standard 
deviation should really be around 0.005. For scoping calculations, perhaps a larger value 
is better, at least until debugging of other aspects of the code can be accomplished. 

7.6.2 Problems with High Resistivity Plasmas 

Two of the most serious problems occur for simulations of high resistivity plasmas 
(“res”~ 1 0 or higher). For instance, in a highly resistive plasma, when there is virtually no 
shock, only a very weak linear wave motion at the beginning of a simulation, then the 
MHD transport algorithm time step is dominated by the magnetic diffusion time dtdif^dx , 
where g=l/res, is the gas conductivity [cm 2 .s] in rationalized Gaussian units. So for high 
resistivity the diffusion time is small, and if diffusion dominates fluid advection (as will be 
the case early on before high speed shocks are fonned) then this time step can be 
significantly smaller than the time step governing advection that satisfies the CFL 
condition. Hence a simulation can advance a large number of time steps yet only a very 
few seconds, maybe only a fraction of a microsecond. (This will often lead to the problem 
of total cell collapse during coarse grid generation mentioned above.) 
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A more serious problem can occur when the mesh size is large — perhaps because one is 
trying to increase the allowable time step governed by time dt^f-qdx' . Notice that when 
the magnetic Reynolds number is large (strong R-field diffusion/weak advection) then dtdff 
must be used to set the simulation time scale for numerical stability. If one does not use 
dtdiff then negative pressures will result from diffusion of the strong magnetic field at the 
boundary of the grid that simulates the shock generator. This may be avoided in future 
updates of the MHD1DSHOCK code by smearing out the pulse current shock source over 
more grid cells. But one can still encounter negative pressures when the shock generating 
magnetic field is turned off. The magnetic field should evolve self-consistently when the 
boundary conditions imposed change from the strong applied shock-driving field to 
continuous boundary conditions. But if the grid spacing is too small then when the 
boundary conditions are changed in this manner, at a single time step, it can cause the 
magnetic field to diffuse too fast for the simulation time step (which is now being governed 
by the advection from the fast shock waves that have been produced). One needs either a 
non-global time step that differs for different regions in the gas depending upon the local 
dominant transport process, or one must protect against negative pressures resulting by 
increasing the coarseness of the grid, making dx larger but for the same shock tube length, 
by reducing the number of grid points, alternatively one could place a multiplicative scale 
factor in dtdff-qdx , such as dtdiff) ■ 1 Qdx , in order to reduce the simulation time step 
artificially. This factor is tunable and is set through “diflf tscale” in the namelist file 
‘mhdini.dat’, so the diffusion time step is defined as dt2=difif_tscale*dx*dx/res, in the 
subroutine SETDT. The actual time tep used in the solver routine STEPON is the 
minimum of this dt2 and “dtl” obtained from the CFL condition on the wave velocities. 
Either fix might avoid negative pressures, but this is not a universal panacea for this 
problem. 

In practice the negative pressures damp out and disappear after a few time steps once the 
transverse applied field at the boundaries is allowed to relax. A problem is that there is a 
serious effect on the fluid variables: One can get gas densities nearly an order of 

magnitude greater for the exact same initial conditions and settings whenever the diffusion 
problem is encountered and the pressures go negative near the boundaries. So the integrity 
and validity of the entire simulation becomes suspect. By changing the length or time 
scales to avoid negative densities one gets (for the exact same inputs otherwise) a much 
lower peak density (and lower peak pressures and temperatures as well). For most of the 
case studies presented in Section 8 of this report, inputs and initial/boundary value 
conditions were chosen to avoid negative pressures. This can be done without too much 
expense by running the MHD1DSHOCK code with calls to MCNP turned off. This runs 
the code very quickly (run times typically in minutes on a Pentium desktop) and one can 
inspect for negative pressures. If the pressures are positive globally then the code can be 
rerun with MCNP turned on, which takes many hours to run. 
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A more permanent fix for these shock tube simulation issues will be attempted in future 
code versions. One solution for example might be to use a smoother shock generation 
transverse magnetic field, one that does not get switched from applied field at the fluid 
boundary to continuous boundary condition at a single time step, but rather a shock driver 
that relaxes over a few time steps. This would be a more pleasingly physical resolution 
after all. However, it does not mean that the essential difficulty with the numerical 
algorithm should be ignored. One would like to avoid negative pressures irrespective of 
how the boundary conditions are applied, as a matter of principle and simulation 
programming hygiene! 

Although not comprehensive, here just a few of pertinent issues that arose, during code 
debugging and simulation testing. The first two deal with getting the simulation to run 
more efficiently and avoiding nonphysical results like negative pressures. The last topic 
here deals with engineering matters on how to simulate greater fission energy release. 

In summary, if you set the input variable “res” to a large value (cm /s) in the file 
'mhdini.daf then you might expect magnetic field diffusion to dominate over fluid 
advection initially in a simulation before strong shocks start to form. This gives a different 
characteristic time dti-O.Sdx !g, that governs numerical stability in the presence of 
diffusion in magnetic field transport. If negative pressures result when g is large therefore 
one should do either or one of at least four things: 

(1) Increase the shock tube length sx, to increase dx, keeping the number of grid cells 
fixed (making a coarser mesh). 

(2) Decrease the number of grid points mx, (again to increase dx by making the mesh 
coarser). 

(Either of these will give larger diffusion times, making it more compatible with the default 
dt\ from the CFL condition. However these are not guaranteed to fix the negative 
pressures; for a fine resolution, one needs some other fix.) 

(3) If you want a fine mesh then you can artificially decrease the diffusion characteristic 
time, say by setting dt 2 =0.1*dx*dx/res, in the source code ‘subroutine setdt’, or similarly 
adapting the scale for diffusion. Because the simulation time step for the leapfrog algorithm 
is set to dt=vcim ( dtx , dtl), where dt\ is the CFL time step from the wave velocities, this 
should ensures a time step that avoids unphysical pressure perturbations from strong 
diffusion. To avoid negative pressure entirely for high resistivity trials one really should 
add extra code to adapt the time step whenever negative pressures are seen, so that the code 
recalculates the transported fluid variables with a smaller time step. (This should be done 
be rewinding the simulation, so that the previous time step variables are not over-written 
before this inspection of the pressures is made). Note that flux-correction should be able to 
avoid negative densities, so generally you should only be worried about the pressure as far 
as problematic diffusion scales become evident. 
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(4) Finally, one can reduce the magnetic Reynold’s number by decreasing the intensity of 
the applied tranverse magnetic field by decreasing the input “byi”. However, if the applied 
shock-drive field is too low then now shocks will form and the purpose of the simulation is 
defeated! So one must in practice run the MHD1DSHOCK code with MCNP turned off to 
perfonn rapid testing of initial conditions to see what sort of magnetic field and pulse time 
is required for adequate looking shock build-up intensities, making sure that the simulation 
runs long enough to capture the shock collision. Then as long as no unphysical negative 
pressures are encountered one can run the script for the coupled codes with MCNP turned 
on. 


7.6.3 Low Fission Power Output 

To boost the fission power production one can either increase the ambient neutron 
population available for multiplication in the pulse (most likely by adding an extraneous 
neutron source, set by the input “source”, which is interpreted by MCNP as ambient 
neutrons per second), or increase the size (physical dimensions of shock tube radius) of the 
reactor, or increase the uranium enrichment or add plutonium (or even add very small 
quantities of americium) to the fuel, which have the effect of increasing k e jf and also the 
neutron multiplication factor in the pulse. 

(1) To increase uranium enrichment, edit the file ‘mcnpfix’, go to the line, 

ml 92235.61c 0.1 92238.61c 0.9 9019.60c 4.0 $default core material 

and change the fraction after 92235.61c to 0.95 say for 95% enriched fuel, then you must 
change the fraction after 92238.61c accordingly to 0.05. 

(2) Increasing the reactor size. To increase the shock tube radius, edit the file 'mcnpfix', 
at the lines, 


begin 

surfaces 

1 

cx 

50 

2 

px 

-50 

3 

px 

50 

10 

cx 

100 

11 

px 

-140 

12 

px 

140 


end surfaces 


that go into writing the MCNP surface card, “1 ex 50” means a radius of 50cm. So to 
model a 160 cm diameter reactor change “1 ex 50” to “1 ex 80”. But then you must also 
give the reflector wrapped around the shock tube more room, so add the difference to the 
number at surface 10, so change “10 ex 100” to “10 ex 130” for the 160cm diameter 
change. 
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(3) Increasing the neutron source. We found that the shock-driven simulation requires a 
large background neutron population for multiplication in the pulse in order to obtain a 
useful amount of fission power. In a real shocktube reactor, the neutron population 
available for multiplication in the pulse will be a result of delayed neutrons from decay of a 
population of delayed neutron precursors generated and sustained by earlier pulses. 
Extraneous neutron sources may also be added to increase this population available for 
multiplication in the pulse. So, to model a background neutron source the parameter 
“source” is specified in the ‘mhdini.dat’ file. For example, set, 


source=l . 0E+10 


in the namelist file ‘mhdini.dat’ to model a background neutron source giving isotropic 
neutron source strength of 10 10 neutrons/sec. If the peak gas density is about 0.1 g/cnr 
then a strong neutron source is recommended to boost fission power production, say: 
source=1.0E+15. 

The simulation records, as a function of time, the system (effective) neutron multiplication 
factor, k e jf. MCNP4C detennines k e ff by Monte Carlo numerical trials, so it is detennined 
statistically from tracking many neutrons through the system. Now, the point of this is that 
it is possible for k ^ to increase well above unity when the fissile fuel density becomes very 
large and a lot of mass is compressed into the reactor core region. This can at sufficiently 
high levels lead to “run-away reactor” and a supercritical chain reaction, which is the 
purpose of the shocktube reactor design since fission power is generated in short bursts 
through the multiplication of existing neutron sources in the pulse. There is no danger 
involved in this scenario since first, it is difficult to achieve the gas densities necessary for 
the super-critical state but also because there is a limit on the total amount of fissile gas 
available in the system. However, more importantly are all the negative reactivity feedback 
effects that will lower k e jf eventually bringing it to a subcritical state before the next pulse. 
These effects are largely derived from the increase in temperature of the gas leading to 
Doppler broadening affecting resonance escape probability but more importantly 
decreasing the density of the gas in the pulse. One potential simulation problem can result 
if the ambient neutron source strength in the shocktube available for multiplication in the 
pulse is set too high. In this case, the multiplication in the pulse can cause the neutron 
population to grow to very large numbers leading to floating-point number overflows in the 
neutron population scale factor “noft” which will crash the simulation. In which case, the 
same input parameters might be selected again, but with lower neuron source strength. 
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7.7 


Shocktube Modeling in MCNP 


7.7.1 Purpose of the Model 

At each step during the propagation of the shockwaves and throughout their interaction, the 
fissile gas and neutron population are contributing to alter the intensive properties of the 
medium through nuclear interactions. At a minimum, it can be understood that these 
interactions affect temperature and molar concentration of various species including fissile 
material, fission products, ionized species, etc. — namely through the fissioning of the 
fissile gas and the ionizations resulting from neutron and fission product interactions with 
the medium. Magneto-gas dynamic modeling of the shockwave propagation and 
interaction is a challenging problem in itself. The additional complexities associated with 
modeling the simultaneous nuclear interactions only make the problem more difficult. 
Therefore, it was decided to develop a series of computer codes to integrate the magneto- 
gas dynamic computer code, MHD1DSHOCK, to model the shockwave propagation and 
interaction and the general Monte Carlo nuclear transport code MCNP5 to model the 
nuclear interactions. These codes are run in series over a specified time step (and later by a 
more physical criterion 14 ) to update the magneto-gas dynamic and nuclear parameters 
before progressing to the next time step (see Figure 1). There is no reason why the nuclear 
code cannot be run after each single time step of the MHD code, but for practical reasons 
this would result in a very long simulation run time, and is probably unnecessary because 
the time scale of the MHD code time step is likely to be much smaller (set by the CFL 
condition) than the typical neutron lifetime. The supporting theory and modeling 
implemented in each code is described below as well as the pre-/post-processor codes that 
facilitate the communications between MHD1DSHOCK and MCNP at each time step. 


14 The MHD code can halt to call MCNP after a fixed time step, but other more physical criteria 
have been implemented as well, so that MCNP will update the fission power and neutron flux 
based upon one or more criteria. Three proposed criteria (which could be used simultaneously) 
are (1) when the maximum velocity in the gas would transport a particle by more than a 
neutron mean free path (this is currently not implemented), (2) when the gas density in any 
cell changes by a significant percentage sufficient to warrant updating the nuclear variables 
(this has been implemented), (3) whenever the elapsed time since the previous call to MCNP 
approaches a few times the neutron lifetime (this has been implemented too). If any of these 
criteria are true then the MHD code will halt to let MCNP run. 
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Shock Tube Model-Nuclear Code 

The shocktube is modeled as a cylinder parallel to x, y, or z-axis with its center point fixed 
at the coordinates (0,0,0). Any differing coordinate system utilized by the 
MHD1DSHOCK code is accommodated by the intermediate code MHD2MCNP that 
prepares the MCNP input based on the MHD1DSHOCK output (for example the MHD 
coordinate system origin is indeed at the left boundary rather than the shock tube center). A 
number of fixed MCNP input cards are supplied in the input file, MCNPFIX, which 
provide the general setup of the problem and are located here for convenience since small 
routine changes and additions can be made without modifying the MHD2MCNP processor 
code. Further, the MCNPFIX file contains all the infonnation on the geometry and 
materials external to the shocktube including the end regions and the radial reflector and 
materials including coils, magnets, and any other structural components. 

The interior of the shocktube is subdivided into cells axially according to the cell structure 
supplied by the MHD1DSHOCK code. Attempts are made by MHD2MCNP to collapse 
the fine cell structure of the MHD1DSHOCK code into a smaller set of cells for use in the 
MCNP calculation (see Figure 5). The criterion for collapsing cells is based on the relative 
change in density from one cell to the next. Specifically, if the maximum relative 
difference between two cells is less than some user specified value, then the two cells are 
marked for collapse according to the following criterion, 


| RhOi ~ Rho i+j 
min (RfiOj - Rho i+ j ) 


< dblMaxRhoDiff 


where i is the left-most cell for collapse and j is a successive cell for comparison. 

Comparison is then made with the next adjacent cell comparing to the original, left most 
cell in this series of collapsing until the maximum relative difference between the current 
cell and the left most cell exceeds dblMaxRhoDiff. Then a new series of comparisons are 
begun starting with the cell with relative density that exceeded dblMaxRhoDiff in 
comparison with the original left-most cell in the series collapsed. So for example if 
dblMaxRhoDiff is 1 .0, cell collapse will continue until the density of a successive cell 
differs from the starting cell by double or half and excluding that final cell where the 
density was greater. 


29 



Reflector and Auxiliary Structures/Materials 



Reflector and Auxiliary Structures/Materials 


(a) 


Reflector and Auxiliary Structures/Materials 



0 >) 

Figure 5. Shocktube models used in the combined simulation codes. 

a) Illustration of the shocktube model output by the MHD1DSHOCK code with many cells, b) Illustration of the 
shocktube model after cell collapse by MHD2MCNP where many of the cells with nearly the same 
density/pressure are collapsed to simplify the MCNP model 

By decreasing the parameter dblMaxRhoDiff one then gets finer and finer coarse meshes 
until eventually the entire MHD grid resolution is resolved, increasing dblMaxRhoDiff 
produces coarser meshes until eventually the routine collapses the entire domain to a single 
cell with density equal to the average over the MHD grid. In this manner an adaptive 
mesh/cell spacing scheme is enabled to capture the fine details of the shockwave profile 
and also reduce the computational burden of tracking many cells in MCNP. 

In collapsing successive cells, the macro cell created has the volumetric average density of 
the micro cells of which it is comprised. The new cell densities along with their cell 
number and location of the right bounding surface are written to a file MHDSIN.DAT for 
later use by the post-processor code MCNP2MHD and the MHD1DSHOCK computer 
code itself for the next time-step. Also, the cell volume and the area of the left and right 
bounding surfaces are calculated and written on VOL and AREA cards added to the MCNP 
input file specification and are required in calculating flux and fission heating tallies for the 
cells and surfaces comprising the shocktube core. 
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Calculation of k eff , Neutron Flux and Cell Heating 

The neutron multiplication factor, k e jf, as well as a number of tallies are made by MCNP to 
predict the neutron flux (per cnr-sec) and heating due to fission energy release (MeV/g or 
MeV/cm 3 ) in each cell. MCNP executes the problem defined by the MHD2MCNP pre- 
processor code and runs for a specified number of histories. The post-processor code 
MCNP2MHD examines the output and detennines if the problem has converged, that is to 
say it has achieved the user specified accuracy for k e ff and the several tallies for flux and 
fission heating. If the maximum allowable standard deviation is exceeded for C,// or any of 
the cell tallies, then the problem is continued for an additional user specified increment set 
from the problem specification file, MCNPPROB, and continues in this fashion until k e f 
and all tallies meet the required accuracy. 

Flux is estimated by a track-length tally (type F4) for each cell defined in the MCNP input 
file. Because the adaptive mesh/cell spacing scheme can produce very thin cells in trying 
to capture the potentially sharp shockwave profile, a surface estimate of the flux (type F2) 
is also tallied at the left and right boundaries of each cell. For very thin cells this should 
provide a better estimate of the flux and reduce the number of Continue MCNP runs that 
must be executed in order to obtain the specified user accuracy. 

Cell heating due to fissioning is estimated using a type F7 tally which measures the energy 
deposited from fission in units of MeV/g for a cell. Alternately, volumetric heating in units 
MeV/cm can be selected as an option on the problem specification file, MCNPPROB. 
Gamma rays from fission are included in this tally by counting them as being deposited 
locally at the site of fission. Additionally, a modified surface flux tally (type F2) is used to 
obtain fission heating estimates for very thin cells as discussed above for flux estimation. 
Here a tally multiplier is specified to account for the correct material densities, fission 
cross-section MCNP reaction number, and energy from fission MCNP reaction number. 

Post-processing of MCNP Output and Preparation for Next 
Time-Step 

The MCNP2MHD post-processor code examines the MCNP output after the initial and 
each subsequent Continue MCNP run (if any) to check the accuracy of the estimates for k e f, 
flux tallies, and fission heating tallies. If the estimated error is within the user-specified 
limits for each, then flux and fission heating tallies are extracted from the MCNP output 
file and appended to the MHDSIN.DAT file initiated by the MHD2MCNP pre-processor 
code. Additionally, the neutron population for this time-step is updated based on the new 
k e jf and prompt neutron lifetime (/„) calculated by MCNP for this time-step and according to 
the following fonnula, 


N(t n ) = N(t n _ l )e (k "~ 1)At/l '’ +^ / ^(e a "“ l)A/// " -1) 

k n ~ 1 

where, n is the n th time-step, At, and t n =t n .i+At; S is the background neutron source strength 
in neutrons/sec. 
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This neutron population is used by the MHD1DSHOCK code to scale the flux and fission 
heating tallies since each of these are normalized to 1 starting neutron. Because of cell 
collapsing in MHD2MCNP, the cell structure of the MHDSIN.DAT file will differ from 
that originally output by the MHD1DSHOCK code in MHDSOUT.DAT and used in the 
subsequent time-step, interpolation on the coarser MCNP grid is required. As a result of 
the interpolation scheme used by the MHD1DSHOCK code, a minimum of six 
interpolation points/cells are required in MHDSIN.DAT. 

One reason for coupling the MHD-gasdynamic simulation code with the MCNP neutron 
transport and fission power tally code is precisely because the neutron flux depends upon 
the multiplication of the neutron population N(Y) in the previous time step, so one needs the 
full simulation history to detennine how the fission power will vary. That is why one 
cannot simply run the MHD shock code and capture frozen instants from the simulation 
and run an MCNP calculation for that instant to determine fission power output, one would 
not know what absolute neutron number density to use as the input in such a scenario. 
These isolated criticality simulations give neutron flux and fission heating relative to one 
source neutron. So there is really no escaping running the codes as a coupled simulation. 
The only alternative would be to build up a massive database of neutron flux and fission 
power as functions of many gas density profiles and histories. This is almost certainly a 
task that would take as long as or longer than simply running the two codes together in 
series with each updating input data for the other as time is advanced. 
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7.7.5 Preliminary Investigations 

A preliminary investigation was made into the methods available to characterize the reactor 
core. To begin with, a simple model of a reactor core with one homogenous region was 
setup in MCNP. This was modeled as a 3 m right-circular cylinder with an additional 40 
cm BeO reflector (see Figure 6 and Figure 7). The fissile gas was assumed to be UF 4 at 5 
MPa partial pressure and 4000 K. The temperature was selected partly based on parsimony, 
since it requires excessive work to incorporate temperature dependent cross sections into 
MCNP. 15 It was first sought to characterize the neutron spectrum in various regions of the 
reactor. The core was subdivided into three regions from the inner to the outer core region 
and the reflector into two regions, inner and outer. To get a detailed spectrum, 482 energy 
groups was setup by thinning the energy grid from the 235 U library used in the calculation. 
New high temperature libraries processed from NJOY99 were used for this evaluation. 
The MCNP kcode was setup to run 100,000 histories in each of 200 cycles for a total to 20 
million histories. Track length estimations of neutron flux were made using tallies for each 
subdivided region. The preliminary calculations showed that the spectrum becomes 
slightly harder moving toward the core center as expected. Detailed examination of 
smaller ranges of neutron energy showed a rather poor characterization of the flux 
particularly near resonances. For this reason the energy grid will be reevaluated to 
redistribute or add points in these regions for future models. 


15 For a full simulation the MHD code will keep track of the true gas temperature that effects 
density and pressure. If required at a later stage of research then variable temperature effects 
(Doppler broadening of resonant for fission and capture) can then be incorporated, though it 
should be noted that this is not a trivial task and requires significant computer memory to be 
permanently allocated for the cross section libraries. 
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Some interesting but unexpected results were obtained when modeling the criticality of the 
shock tube as described in the earlier sections above. A negative reactivity feedback effect 
was noted for two approaching shockwaves. This was modeled for the simulated case of 
two shockwaves of unvarying characteristics as they approach each other from opposite 
ends of the shock tube. That is to say that the strength and profile did not change until they 
interact at the center of the shock tube. It was only then when they interact and collide that 
a positive reactivity is noted and the shock tube core went prompt critical. The reason for 
the negative reactivity associated with approaching waves relates to the effect of the 
reflector region and is analogous to the self-shielding effect noted for the lumping of fuel to 
form a fuel/moderator lattice and can be explained as follows. The rather thin/flat disk 
geometry of the shockwaves is such that the leakage from a shock wave is large in any 
configuration. When the shock waves are far apart, there is a high probability that neutrons 
escaping a shock wave will enter the reflector and slow down in energy before being 
absorbed in one or another shock wave there by increasing the resonance escape 
probability. However, as they approach each other, the probability increases that neutrons 
leaking from one wave will be absorbed in the other wave before slowing down. It is only 
when the waves collide and experience a non-linear increase in density/pressure do we see 
an increase in k e jj since the more compacted, collided waves have a smaller total leakage 
tenn than the two waves separately. This negative reactivity behavior of approaching 
shock waves might be altered by adding a moderating species to the core in the form of a 
low atomic number gas such as CF 4 . However, these and other possible studies are beyond 
the current scope. 8 

Early progress was made with detenninistic codes that are being used to supplement 
current MCNP analysis to reduce uncertainties in the neutron dynamics. This effort has 
proceeded to setup detenninistic modeling of the reactor using the 3D discrete ordinates 
transport code, DOORS3.2. This is to complement some of the Monte Carlo modeling 
perfonned earlier and in the present simulations of Section 9. 

The 3D graph in Figure 8 is of a quarter section of the reactor core that for now is assumed 
to be a 3 m right circular cylinder surrounded by a BeO reflector. The peak in power 
production is at the axial center of the core over an assumed 2 cm width of high-density 
fissile gas constituted from the interacting shock waves. The intent here is only to capture 
what a typical power distribution shape might look like for a generic PMI-GCR shock tube. 
Thus, in particular, the magnitudes indicated in this plot may vary for a given set of reactor 
parameters. The relative values however are meaningful and could be used to gauge the 
power distribution for the ID simulation results to be presented below in Section 8. 


34 





Figure 7. Top view of the MCNP model of the simple core and reflector. 
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A multigroup library was generated for these calculations using ENDF/B-VI neutron data 
and the NJOY99.81 computer code. This work to produce a set of high temperature 
multigroup libraries is progressing. We have noted a couple "bugs" in the code as some of 
the results failed a "sanity check". In particular, NJOY failed to produce realistic 
multigroup values for (v Of) at low energies. We believe this can be attributed to the fact 
that these low energies are represented by two constant vectors at these low energies. To 
provide a temporary fix and produce the attached graph, some hand calculated v Of values 
for the low energy groups was provided. Current and future work will investigate other 
codes/methods to work around this problem as additional work is done on producing these 
libraries. 
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Figure 8. Preliminary results of the fission interaction rate. 
The shock tube radius is 80cm, and the moderator thickness 40cm. 


The intent is not to show any hard data here, it is only to capture what a typical power 
distribution shape might look like for a generic PMI-GCR shock tube. Thus in particular 
the magnitudes indicated in this plot should be ignored. The relative values however are 
meaningful and could be used to gauge the power distribution for the ID simulation results 
to be presented below in Section 8. 
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8 RESULTS OF THE MHD-MCNP SIMULATION 

The results are presented below as three case studies. They have been catalogued by (i) the 
assumed electrical conductivity c !{q 4n) of the gas (low or high), (ii) the type of shock 
generation (“compressor mode” or equivalently “MPD mode” referring to a fixed step 
current at the boundary that is turned off as soon as the shocks collide in the center of the 
shock tube, and “pulsed mode” referring to a step current at the boundary that is turned off 
after a fixed number of time steps), and (iii) various sub-cases depending upon things like 
shock tube length and/or radius variations, fuel enrichment, initial pressure an density, and 
so on. In practice the distinction between shock generation methods “compressor” and 
“pulsed” becomes blurred and indistinguishable when the pulse time is set fairly large and 
the current density at the boundary is reduced for in that case a pulse becomes almost 
indistinguishable from a short duration MPD compressor mode. The term “MPD 
compressor mode” is used because for smaller current densities and longer pulse duration 
times the shock generation arises mainly from nonlinear wave steepening as described 
above and in reference [16 l For short pulse times and large current densities the shock 
generation can be said to arise more from steep magnetic pressure gradient at the 
boundaries giving rise to shocks from the feedback of magnetic pressure into the gas 
through classical magnetogasdynamic coupling. There is however a regime where a true 
pulsed mode could be identified as distinct from a compressor mode by examining the time 
scales for collisional, relaxation and ionization processes, such an analysis should reveal 
when the transition from shocked induced gas to induced waves occurs. This report did not 
look closely at these issues however, so the pulse mode/compressor mode distinction used 
here to characterize the different studies is somewhat arbitrary. 

Before looking at the MHD shock simulations it is instructive to see what happens when 
two gasdynamic shocks collide. 

The main thing to note is that the density sums more or less linearly, as in wave 
superposition. Some nonlinear effects are evident but viscous dissipation and thennal 
diffusion spread the shocks out sufficiently for most nonlinear effects to be washed out. In 
general this is not desirable from a nuclear design point of view because one wants the gas 
density (number density of fissile atoms) to increase by orders of magnitude to ignite 
prompt fission supercriticality. 


[ 16 ] C.C. Wu, “Formation, Structure, and Stability of MHD Intermediate Shocks”, J. Geophysical 
Research, 95 , No.A6, 8149-81475, 1990. 
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It appears then that, from the hydrodynamic shock results, the only way to design a useful 
pulsed reactor would appear to have a strong ambient neutron source placed near the 
reactor to enhance fission power production in the pulse. We will see that this is the same 
for fast MHD shocks with one notable difference, the fast MHD shocks will sum 
nonlinearly when they collide, but most of the energy goes into mechanical pressure in the 
gas, and the density stays fairly low, although it does increase by a few orders of 
magnitude, just not enough, in the absence of a strong background neutron source 
available for multiplication in the pulse, for massive (mega-joule) fission power output. 
This is a somewhat technical distinction that nuclear engineers will be most familiar with, 
which is to say that subcritical systems are generally designed to achieve high thermal 
power output by using an external neutron source to boost the total number of fission 
events, not to necessarily increase the neutron multiplication factor k e fj, which may remain 
low or at some arbitrary value. The gas does achieve a supercritical state however in the 
PMI-GCR shock tube simulations, so the basic viability of the concept is proven though not 
in the way one might have envisaged (one would have expected an intense pulse of power 
rather than a monotonic increase roughly following the shape of the gas density, see the 
results below in Sections 8.1 to 8.4.). The remaining question concerns efficiency and 
practical engineering constraints. 
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Figure 9. Collision of two hydrodynamic shocks in UF 4 gas. 


Next, some case studies of various parameter settings and their effect on fission power 
production and shock collision are presented. The results to follow are obtained by running 
the MHD1DSHOCK code back-to-back with MCNP4C as described in Section 7.4. The 
shocks that will now be investigated for the remainder of this report are fast MHD shock 
waves that carry a strong transverse magnetic field, thus, they are so-called “switch-on” 
shocks. 

For each case study presented, a table is provided showing first the geometry and fuel 
enrichment, and then the initial conditions and boundary conditions driving the flow 
(represented by the applied transverse B-field at the boundary and it’s pulse time). 
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8.1 


Case 1-Low Resistivity Plasma, Compressor Mode 


Case 1-a Results, Low Resistivity, Strong Source 

The first case here had a 100 cm shock tube length, core radius=80 cm, BeO reflector 
radius=130 cm, reflector thickness=40 cm. Fuel enrichment=95% 235-U. The resistivity is 
an almost ideal low 10 111 [cnr.s '] (in rationalized Gaussian ‘cgs’units) throughout. This 
type of core is a thick pancake shape (diameter exceeds length) and the fission power is 
large throughout the length due partly to the BeO moderator extending the length of the 
fuel and beyond. This is perhaps not the best configuration for a well-defined pulse of 
fission power at the local shock collision region, but with the shock driver in MPD 
compressor mode and with a strong neutron source, fission power densities of 
~7xl0“ MeV.cm or — 1 . 14 J.m were obtained after the shock collision. Neutron flux 
levels also climbed to 1 (f 4 n.cm 2 .s '. Although these numbers are impressive one must 
remember that this case used an irregular core geometry that might not be conducive to 
realistic shock tube operation. The large diameter to radius has never been used in an 
electromagnetic shock tube before. 


Table 1 Input parameters for Case 1-a. 


Core length 

Core radius 

BeO reflector 
thickness 

Reflector 

radius 

Grid spacing 

Fuel enrichment 

100 cm 

80 cm 

40 cm 

130 cm 

0.167 cm 

95% 


Po 

Po 

T 0 

By (right) 

Pulse time [ms] 

Source 

[dynne/cm 2 ] 

|g/cm'| 

[°K] 

[Gauss] 

(approx.) 

[n per s] 

10 7 

0.024 

1567 

+708 

1.1 

10 15 


40 



Frames taken from the time evolution of this simulation are shown over the page for gas 
density, fission power density and neutron flux. The results seem somewhat 
counterintuitive in many respects. In particular the gas density peaks when the shocks 
collide but then only very weak outgoing density waves are reflected. The pressure and 
temperature for the outgoing waves are not extremely high — the temperature reaches a 
peak of 3500°K, but this is without the massive heat generation from the fission power 
source added to the MHD equations. So in actuality one would expect extreme fission 
heating of the gas and a very strong outgoing blast wave that would almost completely 
ionize the UF4 gas mixture. This is perfect for conversion to electrical power with a MFCG 
or other inductive circuit. Nevertheless, even with the fission heating source term 
uncoupled, the behavior of the pure MHD shock collision is curious and interesting. One 
thing should be noted, which is that the way the shock driver is set up the shock tube must 
have a make-up supply of gas fuel in order for fluid to be compressed and built up in the 
core. In other words, the MPD compressor river sucks in (or pushes in) a huge amount of 
gas during the simulation (i.e. mass is certainly not conserved within the simulation control 
volume, or more exactly, total control volume fluid mass is not a constant during the 
simulation because of the large (real) mass flux across the boundaries). There is a huge 
pressure built up that does not allow the gas to escape rapidly. It may be possible to 
somehow engineer the system to recover a more classical shock collision and reflection, 
but without the fission heating coupled to the MHD equations it is premature to speculate 
upon what might be the most favorable configuration. 

The fission power density follows the spike in gas density only weakly and is otherwise 
fairly uniform throughout the core. This can be understood because of the large effect of 
the reflector. The mean free path of the fast neutrons bom in the reactor is long because of 
the relatively low density in the core compared with condensed matter. Many such 
neutrons that leak into the reflector slow down in energy before being scattered back to the 
core. These thennal or nearly thennal neutrons have much shorter mean free paths and 
experience absorption or fission largely in areas near the walls of the shock tube. Since the 
current simulation is one-dimensional this radial dependency is averaged out. What is 
observed is the larger amount of fissions occurring in the shockwave because of its greater 
density of fissile material and neutron flux but also to a lesser degree in adjacent cells since 
these are in proximity to the same thermal neutrons emerging from the reflector. This 
explains the somewhat broadened peak in fission heating. 

Also noticeable is the fact that the fission power and neutron flux continually increase as 
time advances following the shock interaction. This again is explained by the huge amount 
of fluid that is transported into the tube and compressed by the boundary current and 
transverse magnetic field pressure. 
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The fission power density axial (x-coordinate) spatial variation follows the gas density 
roughly, though it is smeared out somewhat. The neutron flux appears quite flat however. 
This is quite natural. What is plotted in the case of the neutron flux below is the integrated 
neutron flux over all neutron energies. Referring back to Figure 8 in the previous section, 
one can see the axial and radial dependence of the fission power density for a mock-up of a 
shock collision profile. In that figure one can readily reinterpret the vertical axis as the 
thermal neutron flux instead, because the power density correlates most strongly with the 
distribution of thennalized neutrons. The GCR is not a pure fast or slow reactor in this 
sense, the power is produced by both fast and slow neutrons, but clearly the thermal 
neutrons dominate in fission power production. In the integrated neutron flux plots 
recorded below one sees the effect of both fast and slow neutrons on the total neutron flux, 
and this is clearly almost dead flat across the shock tube, because fast neutrons are 
generated throughout the gas fuel and are reflected back into the gas unifonnly as well. 
The reason the power density doesn’t appear so flat is because fission events, on average, 
occur preferentially near the reflector, i.e. in the gas near the walls where thennal neutrons 
reenter, and also near the highest density regions near the shock waves. The mean free path 
of the fast neutrons in the gas fuel is so large that fast neutrons don’t really see the variation 
in gas density. This explains why the integrated neutron flux is so flat. 

The UF4 gas pressure is plotted after the series of neutron flux plots for the same time 
frames as the gas density and fission power density. One can see the rise in gas pressure 
that follows the shocks. However it is clearly not the mechanical pressure alone that 
prevents most of the mass escaping with the outgoing reflected shock fronts. 

Following the series of pressure plots, one can see the same time frame variation in the 
transverse magnetic field that travels with the shock waves. The By component increases to 
a large magnitude and also spikes as the shocks collide. After collision however the rise in 
the By is sustained, but more importantly it is propagated along with the outgoing reflected 
shock fronts that are weaker than the incoming shock fronts. Nevertheless, this strong 
magnetic field exerts considerable pressure on the ionized gas; this explains why the pocket 
of high-density gas remains in the center of the shock tube. 

The gas temperature is also increased by the shock collision; this will raise the electrical 
conductivity of the gas and help to create the same high magnetic pressure. For a 
realistically high resistivity UF4 plasma (with little or no seed material to boos 
conductivity) there is good reason to expect that even for realistic partial plasma or weakly 
ionized gas models the same effect will occur. 
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Case 1-a time frames: Low Resistivity Plasma, Strong Source, BeO End Wall. 
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Case 1-a Fission Power Density, for the same time frames. 
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Case 1-a Integrated Neutron Flux, for the same time frames. 
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Case 1-a Transverse Magnetic Field (b y = B v [Gauss]/ V4?r ), for the same time frames. 
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Case 1-a UF4 Gas Temperature, for the same time frames. 
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8.2 Case 2-High Resistivity Plasma, Compressor Mode 

In this case much higher resistivity plasma was chosen, this would simulate a weakly 
ionized gas. Although not very realistic (because the intense shock waves can almost 
completely ionize a gas at the shock front) this does give one an idea of how a non-ideal 
working fluid behaves. This example is slightly atypical in that the grid cell size was 
coarsened considerably to 1 .Ocm in order to avoid negative pressures forming at the grid 
boundary. A fairly strong applied magnetic field was required to get sizeable densities at 
the shock collision. 

Case 2-b Results, High Resistivity, Strong Source 

Input parameters used are summarized in Table 2. Results are plotted for a few time frames 
over the page. The resistivity is 10 5 [cm 2 . s' 1 ] throughout. 


Table 2. Input parameters for Case 2-b. 


Core length 

Core 

radius 

BeO reflector 
thickness 

Reflector 

radius 

Grid spacing 

Fuel 

enrichment 

100 cm 

50 cm 

40 cm 

100 cm 

0.167 cm 

95% 


Po 

Po 

T 0 

By (right) 

Pulse time [ms] 

Source 

[dynne/cm 2 ] 

[g/cm 3 ] 

[°K] 

[Gauss] 

(approx.) 

[n per s] 

to 7 

0.024 

1567 

+35450 

0.10 

10 15 


The most noticeable difference between this case and the low resistivity plasma (Case 1-a) 
is the speed of the shock. This is produced by the larger applied magnetic field (note that 
bv0r= 10000 “ gaus s” in the figures actually corresponds to a real applied magnetic intensity 
of B v = by 0 44 k [Gauss] =35450 Gauss). The fission power density is lower, partly 
because the neutron population did not have sufficient time to build-up to significant 
numbers. However, there was still a large fission power release, though probably not 
enough for a big power gain. 
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Case 2-b time frames: High Resistivity Plasma, Strong Source, BeO End Wall. 
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As with Case 1-a the density stays large in the middle of the shock tube. The other 
difference in this case is that the reactor is smaller; it is a right cylinder of diameter 100 cm. 
So there is less fuel to begin with in Case 2-b here. The neutron source strength was the 
same and the moderator dimensions the same as Case 1-a, so Case 2-b suffers from smaller 
reactor size and overly rapid shock collision. 

The moderator material at the ends of the shock tube gives a pronounced higher fission 
power density near the end boundaries, although the power rise still follows the variation in 
density fairly closely. For this simulation the fission power continues to rise throughout the 
reactor as time goes by and rises to about 4x10 MeV.cm after 1.60 ms, which is 
considerably higher than the power released at the moment of greatest density. Moreover 
the trend indicates that the power would continue to climb, so there appears to be almost 
enough mass now in the shock tube to generate a run- away surge as happened in Case 1-a. 
This is an effect of the quasi-steady density pile-up in the shock tube that has forced so 
much fuel mass into the moderated region. This points to the need to engineer some kind 
of fuel dispersal after the shock collision (only if the fission heat source doesn’t achieve a 
blow-away of all the mass, which one would probably expect!). 
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8.3 


Case 3-Moderate Resistivity Plasma , Compressor Mode 


Case 3-a Results, High Resistivity, Long Pulse Duration, 
Strong Source 

For this case a time in seconds was used to control the applied field drive time, thus in the 
input table the pulse time is marked as exactly 0.1 milli-seconds. Other parameters are 
typical of those used elsewhere, though the source strength was moderate at 10 12 n/sec. 
The fission power is plotted in units of MeV/cm . This shock tube was surrounded by a 
BeO annulus reflector and had 40 cm thick structural material (Zr) as the end walls, so the 
overall length was 280 cm. The resistivity is 1.0 [cnr.s '] throughout. 


Table 3. Input parameters for Case 3-a. 


Core length 

Core 

radius 

BeO reflector 
width 

Reflector 

radius 

Grid spacing 

Fuel 

enrichment 

200 cm 

60 cm 

200 cm 

110 cm 

0.2 cm 

95% 


Po 

Po 

To 

B v (right) 

Pulse time [ms] 

Source 

[dynne/cm 2 ] 

[g/cm 3 ] 

[°K] 

[Gauss] 

(exact) 

[n per s] 

10 7 

0.024 

1567 

+7090 

0.1 

10 12 


In this simulation the fission power again gradually increases monotonically throughout the 
reactor shock tube core, it climbs reasonably steadily to -3x10 MeV.cm after 2.1 ms, 
and shows no sign of decreasing. 

Accompanying the plots of the gas density is a time history of the neutron scale factor 
(number of neutrons) used by the nuclear code to nonnalize statistics. Over the page are 
the corresponding plots of the fission power density. 
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Neutron Population vs Time 



Figure 10. Neutron population versus time for the long duration pulse (Case 3-a). 


54 


(b) Fission power density. 
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8.4 Case 4-High Resistivity Plasma , Pulsed Mode 

For this study the applied magnetic field representing the shock generating current density 
and ponderomotive force was reduced in duration. 

Case 4-b Results, High Resistivity, Pulsed Mode, Strong 
Source, BeO End walls 

For this case we plot the fission power in units of MeV/g, and place the plots next to the 
corresponding density at the same time. So density is on the left with specific fission 
power on the right. For the input parameters notice the pulse time here of only 0.02 ms. 
The core has a fairly large radius=80 cm to get slightly more fuel mass initially in the core 
to compensate for the lower drive power. This shock tube was entirely surrounded by 
40 cm thick BeO moderator. 


Table 4. Input parameters for Case 4-b. 


Core length 

Core 

radius 

BeO reflector 
thickness 

Reflector 

radius 

Grid spacing 

Fuel 

enrichment 

100 cm 

80 cm 

40 cm 

130 cm 

0.167 cm 

95% 


Po 

Po 

To 

By (right) 

Pulse time [ms] 

Source 

[dynne/cm 2 ] 

[g/cm 3 ] 

[°K] 

[Gauss] 

(approx.) 

[n per s] 

10 7 

0.024 

1567 

+35450 

0.02 

to 15 


The results are not dissimilar to the compressor mode studies. Remember that now the 
specific power is plotted. If the power density was plotted then the results would look very 
much like the history in Cases 1-a and 2-b, only different characteristic times and power 
magnitudes are involved. 

Notice that this means the core (shock tube) is quite pancake shaped, (diameter exceeds 
length). This helps get more mass into the core at high density while limiting dissipation of 
the shock. However it looks as though the shocks do not significantly dissipate, nor would 
they be expected, given the absence of viscous and heat radiation or conduction losses in 
the governing equations. In future simulations however it may prove that a real (viscous 
weakly ionized) gas will produce rapidly attenuated shocks that weaken into linear waves 
and fail to produce high density in the reactor. Such effects might not result in loss of 
power, however, if run in MPD compressor mode, because as the Cases 1 and 2 show, one 
doesn’t necessarily need supercritical densities in confined regions to generate lots of 
power, a smeared out shock or sum of waves might do just as well (ignoring for the present 
time issues of power extraction, that is). 
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Case 4-b time frames, Pulsed Mode, High Resistivity, density and q/ lss plots. 
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Case 4-d Results: High Resistivity, Pulsed Mode, Medium 
Source, Zr End Wall 

Another trial with weaker neutron source strength was simulated. This time with 40 cm 
thick Zr structural walls at the ends of the tube. The shock tube was made fairly large to 
compensate, and the applied magnetic field was increased slightly to about 21300 Gauss, 
but the pulse duration was extremely short, only 5.0 microseconds. For this case we again 
plot the fission power in units of MeV/g, and place the plots next to the corresponding 
density at the same time. 


Table 5. Input parameters for Case 4-b. 


Core length 

Core 

radius 

BeO reflector 
thickness 

Reflector 

radius 

Grid spacing 

Fuel 

enrichment 

200 cm 

80 cm 

50 cm 

130 cm 

0.5 cm 

95% 


Po 

Po 

To 

By (right) 

Pulse time 

Source 

[dynne/cm 2 ] 

[g/cm 3 ] 

[°K] 

[Gauss] 

(approx.) 

[n per s] 

10 7 

0.024 

1567 

+21 300 

5.0 ps 

10 12 


The waves are no longer acting nonlinearly and they are fairly weak. One would expect 
that megagauss fields would be needed to produce strong shock waves with such short 
applied field pulse duration, and this is indeed what these results seem to indicate. Figure 
1 1 shows the neutron population as a function of time. 
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Neutron Population vs Time 


High Pressure UF 4 E-M Shock Tube, 'Pulsed Mode" 



Figure 1 1 . Neutron population versus time for the short pulse trail (Case 4-d). 


The curious thing is that here we see that even for a long time period after shock collision 
the mass of fissile fuel dragged into the shock tube confines can eventually lead to a huge 
multiplication of the ambient background neutron population. Thus, fission energy is 
eventually released, just not as violently, nor perhaps as usefully, as might occur in a classic 
shock collision followed by reflection. 
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Case 4-b time frames: Microsecond Pulsed Mode, High Resistivity, density and q/ lss plots. 
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8.5 Other Case Studies 

A whole series of trials with various input parameters has been under investigation to 
systematically search for an optimum configuration of geometry and shock generator that 
can maximize the peak density and lead to high fission power density during shock 
collision. This is the time at which one wants fission energy to peak so that power can be 
extracted from the explosive effect on the ionized gas. However, to date although most 
trials eventually are found to result in a build-up in neutrons and consequent monotonic 
increase in fission power output, efforts are still on-going in finding either a classical shock 
collision and reflection regime with sudden fission energy burst, or to find a regime where 
the fission energy at least spikes momentarily as the shocks collide. 

There are so many parameters controlling the PMI-GCR shock wave driven system that it 
is prudent to avoid speculating upon all the possible avenues open for improving upon the 
research presented in this report, or this report would become a behemoth. Take just one 
example, if the persistent build-up of fluid in the shock tube is a real effect then it offers 
alternative ways to design a power conversion device from this concept. Perhaps one can 
even start with a very rarified low density, low pressure fluid and still get massive 
quantities of fission energy over slightly longer time intervals, merely by 
electromagnetically compressing the fluid. Such case studies await future endeavors and 
may even constitute a totally new class of competing shock driven reactors. 


9 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 
9.1 Summary of Research Findings 

It has been seen that pure hydroshocks in UF4 dissipate too rapidly for nonlinear increases 
in the gas density required for effective fission power production. With fast MHD shocks 
however it was shown that waves can be formed that steepen into shocks that sum highly 
nonlinearly, but with most energy locked up in mechanical pressure, although the plasma 
density does increase by orders of magnitude as well, only not as greatly as the hydrostatic 
pressure. Yet this may not be enough to achieve a highly supercritical core that would be 
necessary for massive (~ MJ/cm 3 ) fission energy output. The gas does achieve a 
supercritical state however, so the basic viability of the concept is proven. There is certain 
potential for extraction of net power from the PMI-GCR system. The remaining question 
concerns efficiency and practical engineering constraints. In further simulations, release of 
hundreds of thousands of Joules over a few microseconds is possible through the 
introduction of a strong neutron source for significant fission power through multiplication 
in the pulse due to the shocktube prompt critical state afforded by the nonlinear increase in 
density. 
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The recommended gross design parameters that emerge form the research so far, for a 
viable PMI-GCR power supply, are listed in Table 6. 


Table 6. Recommended design parameters for 
an ambient-subcritical electromagnetic shock- 
driven PMI-GCR. 


Core radius 

>80cm 

Shock tube length 

200cm 

Moderator thickness 

50cm 

235 U enrichment 

95% 

Neutron source 

10 10 — 10 15 n sec -1 

Ambient gas pressure 

1 MPa 

Ambient gas density 

0.024 g em -3 

5-field strength (MPD 
mode) 

0.2 to 3.0 Tesla 

.5-field strength 
(pulsed mode) 

20 to 200 Tesla 


Comparison of results for Case 1 and Case 2 shows that higher resistivity partial plasma 
reduces the effectiveness of fission power production (although only indirectly through the 
shock collision effects). It may be that the gas conductivity used for Cases 2 and 4, 
£•=10 " [cm .s] is too low and not representative of the conductivity that a shock heated 
gas might approach. Undoubtedly more investigation needs to be done to ascertain 
whether or not the reactor parameters recommended here will yield orders of magnitude 
lower or higher fission power production. One can only conclude that the reactor is viable 
and that there is a highly novel and interesting phenomenon that occurs leading to sustained 
accumulation of fuel mass in the center region of the electromagnetic shock tube, this 
critical mass allows sustained growth in the neutron population and hence in the fission 
power yield over many milliseconds. So far simulations have not looked forward far 
enough in time to see if this dynamic effect is sustainable. So the PMI-GCR concept is both 
intriguing and worrisome. It is not yet known whether it can be made efficient and 
practical. That is, one cannot yet say whether the system will perform as a useful power 
supply even though under the high magnetic field compression perfonned in the 
simulations above some massive yields of fission thennal power have been seen. 

If the sustained capture of mass in the shock tube core is a real effect then it bodes 
extremely well for the promise of the PMI-GCR class of pulsed fission reactors. One 
would have to ensure that some mechanism was in place for removing (pumping out) the 
significant concentration of compressed fuel in the reactor shock tube core to allow both 
cycling of the device as well as limiting damage due to the rapid increase in fission power. 
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It is hoped that once the magnetogasdynamic simulation code is modified to include the 
fission power source term and a better plasma kinetics model that some of the real gas 
effects will come into play to limit the compression effect without destroying the nonlinear 
shock interaction that produces the high fluid densities in the first place. 


9.2 Recommendations for Future Research 

The research findings reported above are strong results. They are however weakened by 
the lack of complete coupling of the nuclear and MHD codes. A full simulation of the 
shock generation and collision leading to fission energy release has therefore not been 
accomplished. Careful groundwork has, however, been laid that will make the completion 
of these unfinished tasks close to a fonnality with just one qualification: the need for 
gathering experimental data on high temperature UF4 non-equilibrium partial plasma 
properties poses a nontrivial task. 

There are three urgent tasks that are suggested as serious omissions that need to be 
carefully examined in any future research on this concept. 

1. The fission source terms must be included in the MHD equation solver, and 
preferably 2-dimensional geometry should be modeled with viscous and heat 
conduction tenns. Radiative losses should also be considered as the UF 4 partial 
plasma is likely to become optically very opaque during shock collision, but 
yet highly radiative due to the enormous combined effect of shock and fission 
heating that may raise the gas temperature up to millions of Kelvin. 

2 . UF4 plasma thermoproperties and transport properties need to be modeled with 
more rigour. It is highly unsatisfactory to have the gas conductivity fixed at an 
arbitrary value somewhere between the equilibrium scalar conductivity of UF4 
and gas discharge conductivity. Experimental data is paramount for this effort. 
Experiments to measure UF4 partial plasma conductivity and viscosity at high 
temperature are in urgent need to complement limited theoretical transport 
coefficient calculations. 

3. The pulsed R-field shock generation mechanism was inadequately investigated. 

Although the MPD-compressor mode modeling also made gross approxi- 
mations such as zero current rise time and high ionization, the effectiveness of 
such a shock generation mechanism can be questioned and needs to be 
carefully evaluated before a huge amount of effort is put into costly experi- 
ments that might fail. Although pulsed high R-field shock generation is 
technically more demanding on the input power requirements, and also quite 
hazardous, nevertheless, it yet may prove to be the only way to generate strong 
shocks in UF4 partial plasma. The reason is that the ambient UF4 is 
nonconducting, or at least very weakly ionized (~1% ionization or lower) at 
ambient 2000°K assumed for the reactor, so the entire effect of an MPD shock 
generator must depend upon gas discharge ionization, which is a local effect 
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and so does not ionize the gas between the shocks, only the shock fronts and 
the gas behind the incoming shocks will be ionized. 

These tasks can be worked on separately to some extent, so there is ample reason to 
expect that a future phase of research can successfully address these concerns given only 
that sufficient human work hours are allocated. If the above tasks can be performed then 
simulations should be able to aid in the design of effective and worthwhile “zero power” 
shock tube experiments with non-fissile UF4 (or with a chemically similar working fluid). 
This would be to validate shock generation and collision modeling, not fission power 
production. At a later stage of research, to simulate fission power production 
experimentally, a strong external heat source could be used to add thennal energy to the 
shock tube experiment as is found effective in conventional nuclear reactor testing 
programs that do not use actual fuel. 

The analysis of the potential for power conversion has not been accomplished and waits for 
future research efforts. One needs to include the strong fission source terms in the MHD 
equations as mentioned above. Stronger post-collision shocks are expected to result when 
this effect is included, with potential for fairly efficient power conversion using the 
magneto-cumulative flux compression principle. For the power conversion stage at least 
one can say that the gas will be highly conducting due to the ionizing shock waves and 
fission product ionization. This will aid with fission-to-electric power conversion efficien- 
cy. 


* * * 

In summary, strong conclusions can be made on the overall feasibility of the PMI-GCR 
concept. There is no doubt that it can function as a viable power source. In some trials 
power densities of 10 10 MeV.cm are seen at the time of maximum density, and this can 
grow to 10 MeV.cm or even greater if the gas is left unperturbed in the center of the 
shock tube. This is a novel effect that was not predicted and has not been fully evaluated as 
a possible viable new means of fission power conversion. It is not yet fully kn own what 
regime of efficiency the system might operate at, or for what time periods, but the cycle 
time should be on the order of a few milliseconds. A limiting constraint might be how fast 
a capacitor bank can be recharged to start successive shocks, which has not been 
considered in this research. Strong recommendations are also made; especially one can 
single out fissioning partial plasma properties for UF4 as urgently in need of experimental 
investigation. 
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A.l 


ADDITIONAL SUBROUTINES FOR COUPLING THE MHD AND 
MCNP CODES 

This subsection contains a brief description of the subroutines added to the MHD code. 
The script codes for MHD post-processing, MCNP preprocessing, MCNP post-processing, 
and MHD pre-processing will not be described here. The MCNP4C code is documented by 
Los Alamos National Laboratory and so is not described here either. So what remains is a 
brief description of the MHD code structure. 

The main new subroutines added to the basic MHD code are, 

SIM_STATUS: Detennines whether the MHD simulation is continuing from a pause to 
run MCNP or starting fresh, opens data and log files for ‘append’ or ‘new’ access as 
appropriate. 

MHDSAVE (ni,nf, filename): Used to store arrays for MHD1DSHOCK while MCNP is 
running. 

MHDRETRIEVE (ni,nf, filename): Used to read back the array data saved by 

MHDSAVE once the script calls MHD1DSHOCK to resume a simulation. 

MHDOUTPUT (ni,nf,ncols,fname ): Writes data to a file ‘mhdsout.dat’ that can be read 
by the MCNP pre-processor, the data recorded includes cell coordinates and gas density at 
cell nodes, neutron source strength, neutron population from last time MCNP was called, 
and some data for MCNP input cards.. 

READMCNP (ni,nf,fname): Reads in coarse mesh cell coordinates, fission power density 

o 

(=”cqfiss” in the source code) q'" [MeV.cm per neutron] (or specific power c/[MeV/g per 
neutron]), and neutron flux (-’cfluxn” in the source code) [n.cm .s per neutron], and 
interpolates and extrapolates these values onto the fine MHD grid mesh (to arrays named 
“qfiss” and “fluxn” in the source code), scaling them by the neutron population scale factor 
“noft”. 

PLOTOUTPUT: Uses user supplied settings to record all data capturing the state of the 
fluid flow, magnetic field, and source terms. The user can specify one of two formats (e.g. 
Tecplot readable format, or a custom fonnat) and can specify a certain time step interval 
between records by setting the variable “nframes”, and can specify the number of cells of 
the grid to record data for by setting the variable “nxdata” (nxdata<nx is enforced, where 
nx=number of grid cells). 
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MHDSTOP(strstop): Halts the MHD1DSHOCK code with the correct return value to 
alert the script what appropriate action to take depending upon the value of the string 
“strstop”. 


A.i.i New Input Parameters to the MHD Code 

The original hnhdld.for’ source code came with an include file ‘mhdinc.inc’ containing the 
common block variables, and a namelist file ‘mhdin.dat’ for entering the main simulation 
parameters such as the maximum time steps, the Courant number, the gas ratio y (=1.07 for 
UF 4 ), and the plasma resistivity, c [cmVs] (named “res” in the source code). To these files 
an extensive number of additional control variables, parameters and other structures have 
been added, so what was originally a 600 line F77 program called ‘mhdld.for’ has been 
modified and extended into 5900 lines of F77 source code that has undergone three test 
version updates and is now called ‘mhdld_4.for’. It can be compiled using the GNU g77 
FORTRAN to C compiler. 

For a user the most important updates are the new variables and parameters in the new 
(version 4) namelist file ‘mhdini.dat’. This is described next. 

The contents of the input file ‘mhdini.dat’ consisting of 39 lines look at present as follows 
(default settings are somewhat arbitrary for some of these, but most are typical), 

&run 

nend=4000 

nse=10 

nsp=10 

cr=0 . 4 

if luid=0 

gamma=l . 07 

res=l . 0E-10 

nsg=10 

lef tbnd=0 . 0 

sx=100 . 0 

tperiod=l . OE-6 

dtmcnp=l . OE-6 

rhocrit=0 . 026 

ncallmcnp=100 

rhoin=0 .0241 

Pin=l . OE+7 

bxi=l . 0 

byi=0 . 2E+3 

nJzof f=0 

ncurzon=500 

physunits= ' cgs ' 

viewmhd= ' student ' 

nframes=145 

nxdata=200 

strlogf ile= ' mhdld . log ' 
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strmcnpf ile= ' mhdsout . dat ' 

strmhdinf ile= ' mhdsin . dat ' 

strmhdsaved= ' mhdsaved . dat ' 

start file= 'start.ini' 

strmirrorf ile= ' mirror . dat ' 

savedir= ' results/ ' 

ncoremtrl=l 

ncoreimp=l 

source=l . 0E+10 

isf luxlim=F 

isantidif f=F 

ismcnp=T 

isdebug=T 

A brief synopsis these parameters and their functions are given in Appendix Table 7. The 
data types for each parameter are omitted here, but should be obvious. 


Appendix Table 7. List of inputs to the MHD code, listing of namelist file ‘mhdini.dat’. 


nend 

= maximum number of timesteps (iterations) of the MHD solver. 

nse 

= counter for the timestep at which a subset of the array data will be written to the file 
‘mhdld.log’. 

nsp 

= counter for the timestep in which output will be echoed to the screen terminal. 

nsg 

= an additional (presently unused) counter. 

cr 

= Default=0.4, the CFL number to satisfy the Courant-Friedrich-Levy condition. 

i fluid 

= 0, uses ideal gas equation of state to get T(P,rho) 

= 1, then uses user supplied empirical equation of state. 

gamma 

= Gas ratio of specific heats. 

res 

= Electrical resistivity of gas as a partial plasma, units [cmVs” 1 ].. 

leftbnd 

= Default=0.0, axial coordinate (cm) of leftmost fluid cell left hand side boundary. 

sx 

= Shock tube length, [cm]. 

tperiod 

= time period [secs] for calls to write array values to output files for later visualization of 
results. 

dtmcnp 

= time increment [secs] at which to pause MHD program to call MCNP via the batch 
script. 

rhocrit 

= minimum density required in the shock tube before the MHD code bothers invoking 
MCNP, i.e., for MCNP to be called some ro(j)>rhocrit is required. 

ncallmcnp 

= option to specify maximum number of calls to MCNP during the simulation, e.g., if 
ncallmcnp=100 and there are nend=6000 time steps, then using ncallmcnp for control 
will pause to run MCNP every 60 time steps. 

rhoin 

= initial density throughout tube [g/cc]. 

Pin 

= initial shock tube gas pressure [dynne/cm A 2]. 

bxi 

= constant value of i? v /sqrt(4*pi) field. Units B x : [gauss] 

byi 

= boundary value for /i,/sqrt(4*pi) which sets the current J : . by setting the tangential 
components of the magnetic field intensity at the grid guard cells. (Very roughly 
by0(r,l) (-1.0, 1.0) gives / Z =5.0E+12 [statamp/cm 2 ] (l,7E+7 [Amp/m 2 ]), for dx 0.166 
cm, for other dx grid size J : will be different because the exact value of J z depends 
upon dx.) 

nJzoff 

= 0 to turn off boundary current when vx(nx/2) gets perturbed, 
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= 

1 to turn off boundary current after fixed time step. 

ncurzon 


A fixed integer time step after which the current at the boundary is turned off (in effect 
only when nJzoff=l) by resetting the boundary tangential magnetic field in the guard 
cells of the mesh to zero. 

diff tscale 

= 

The fraction by which to decrease the diffusion time scale du. In the code this is, 
dt2=diff tscale*dx*dx/res. 

physunits 


Option to specify whether using Gaussain (physunits=’cgs’) ro SI units 
(physunits=’mks’), eg. one can use this to tell various functions whether to expect 
argument values in SI or cgs units and hence do necessary conversions. For example, 
some of the real gas thermoproperties are sometimes most conveniently entered in SI 
units, but the MHD source code proper employs Gaussian (cgs) units for computer 
integration of the governing equations. 

viewmhd 

= 

type of visualization format to use when writing array data output for visualization. 
Presently must be either 'viewmhd', ‘student’ or 'tecplot' 

nframes 

= 

desired number of time steps (out of a possible total=nend) for which array data will be 
recorded, ie. written to the visualization data files ‘ro.daf , ‘vx.daf , etc...... 

nxdata 

= 

maximum number of points (cell nodes) for which data will be written to 'ro.dat', 
vx.dat', etc... 

strlogfile 

= 

specifies name of a log file that records diagnostic information about the simulation. 
The script does not care about this file. 

strmcnpfile 

= 

Default-mhdsout.dat', the name of the file that the script program expects to see data 
from the MHD code relevant to the MCNP preprocessor script. 

strmhdinfile 

= 

Default-mhdsin.dat', file that the script prepares containing results of MCNP output 
that will be read by the MHD code. 

strmhdsaved 


'mhdsaved.dat', contains all the arrays that the MHD code needs to remember for 
continuing the same simulation after temporarily halting to call MCNP. The script does 
not care about this file. 

strstartfile 

= 

'start.ini', used to evaluate simulation status. 

strmirrorfile 

= 

'mirror.dat', used for diagnosing bugs in retrieving data from MCNP and saved arrays 
from the file named by “strmhdsaved”. 

savedir 

= 

'results/', name of subdirectory to store simulation data output record files to. 

ncoremtrl 

= 

constant used by MCNP preprocessor. 

ncoreimp 

= 

constant used by MCNP preprocessor. 

source 

= 

background neutron source parameter used by MCNP (neutrons/sec). 

isfluxlim 

= 

logical constant , if i sfl nx 1 i m — F, then flux limiting is turned off. 

isantidiff 

= 

logical constant, if isantidiff=F, then antidiffusion step ofFCT is turned off. 

ismcnp 

= 

logical constant, ismcnp=T iff MCNP is to be used. 

isdebug 

= 

If isdebug =F, then some lines of the log file output will be suppressed. 


As one can see, most of the parameters are non-essential for the MHD solver algorithm in 
‘subroutine stepon’, and the majority of the inputs are for controlling how MHD ID- 
SHOCK communicates with MCNP4C. 

The main program loop and main additional subroutines will be listed next. The structure 
of the main program loop (mainly focusing on the program calls and control structures) is 
as follows in Appendix Table 8, for clarity subroutine calls and Fortran keywords are in 
blue, control structures are in red, ordinary assignment statements and so forth are in black, 
pseudocode is indicated in black between # marks, and comments are in plain text prefixed 
by exclamation marks. 
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Appendix Table 8. MHDlD_4.for Main Program Loop Structure. 


Program MHD1D ver4 

issimgo=.false. ! Initially assume that a new simulation case is to be run. 

call input ! Read inputs from namelist file ‘mhdini.dat’. 

call sim_status ! Read the file ‘start.ini’ (named by the variable “strstartfile”) to ascertain whether a 
! simulation is resuming after pausing to run MCNP, or if indeed a new simulation is 
! starting, ‘start.ini’ should contain a single record =”start” or =”continue” as appropriate. 

! A key task of this subroutine is to set the value of “issimgo” correctly. 

if (.not.issimgo) then 

call initial ! Set up initial conditions on the mesh only if starting a new job. 

#set fluxn(j) and qfiss(j) arrays to zero.# 

#initialize variables, e.g. nstep=0# 

end if 

# set integers for array maximums ! Not all subroutines need the guard cells, so the min. and 
and minimums to be passed to var- ! max. index passed depends upon the task. 

ious subroutines, and set integers 

controlling data output preferences.# 

if (issimgo) then 

call readmcnp( ) ! Get MCNP data if resuming a simulation. 

call mhdretrieve( ) ! Get saved MHD array data if resuming a simulation. 

end if 

do while (nstep.le.nend) ! Enter main MHD solver loop. 

call setdt( ) ! Set the max allowable time step using the CFL condition, 

call stepon( ) ! Transport the fluid and field variables one time step. This is the FCT algorithm. 
#Perform some diagnostic checks 
e.g. save global maximum p,P,T# 

call bndry ! Reset the boundary conditions as appropriate for new time step. 

if (nstep.eq.l) cycle 

if (#nstep or time is due for recording a frame of data#) then 

call plotoutput(nxskip) ! Record arrays, but only every nxskip grid point. 

end if 

#Check criteria for calling MCNP to ! This sets the variables “isrhocrit” and “iscallmcnp”. 

update fission power source terms# 
if (iscallmcnp) then 

delta=simtime-prevtime ! This is the time since last call to MCNP, it is required by the 

! MCNP pre -processor. 

call mhdoutput(ni,nf,ncols,strmcnpfile) ! Record data required by MCNP pre -processor, 
call mhdsave(ni,nf,strmhdsaved) ! Save arrays and data required by MHD1DSHOCK to 

! restart from previous state when the script starts it again, 
call mhdstop(‘Pause’) ! Halt execution but return a value=l to the OS that indicates 

! to the script that the simulation should continue after MCNP is done 
! updating the fission source terms. 

end if 

end do 

#Write a summary of simulation results 

to the log file and to the screen terminal.# 

call mhdstop(‘End' 1 ) ! Stop the simulation for good by sending return value-0 to the OS. 

end 
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There are of course numerous other subroutines and user defined functions that are called 
by the main subroutines listed in the table, but it would be pedantic to describe them all 
here. The main critical background subroutine is “Coarse2FineInterp( )” that gets called 
twice by subroutine “readmcnp( )”. Subroutine Coarse2FineInterp( ) takes array indices 
and array data (representing some function on a mesh) as arguments and interpolates data 
defined on a coarse mesh onto a finer mesh. In the MDH1DSHOCK code the coarse mesh 
array is read in from the MCNP results, while the fine mesh is just the MHD grid. Of all 
new procedures in the code this one is perhaps the least well tested and may still contain a 
few bugs. The subroutine “stepon( )” that integrates the MHD equations forward by one 
time step is the other procedure that can be improved, currently it operates best when flux- 
correction and anti-diffusion are turned off, and thus it operates as a leapfrog trapezoidal 
solver with strong artificial diffusion. This will not give unreasonable results for the 
reasonably high viscosity UF 4 gas fuel, but one should bear in mind that the equations 
being solved are not really physically accurate, although they are still qualitatively useful 
for suggesting future experiments and code modifications. 

There are also additional user defined functions that can be used to compute the real gas 
properties as a function of temperature and pressure, however these are generally only valid 
for T=200 to 10 000°K, and for equilibrium UF 4 plasmas, so caution should be exercised in 
using them, and generally it was found that using ideal gas laws supplemented by a fixed 
electrical conductivity yielded acceptable results for the purposes of this study. 
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